一种飞机大尺寸曲面测量点差异性规划方法

2020-06-16毛喆李泷杲徐岩曾琪主逵

北京航空航天大学学报 2020年5期

毛喆，李泷杲，*，徐岩，曾琪，主逵

（1.南京航空航天大学机电学院，南京210016； 2.深圳市劲拓自动化设备股份有限公司，深圳518101）

现代飞机高机动、高气动的产品需求对装配质量优化提出了巨大挑战，数字化测量作为飞机装配中的重要环节，测量质量的优劣性严重影响飞机的实际装配结果。目前高精度、高效率的数字化测量技术正逐步代替传统的检测方式，但测量实施前需基于模型进行测量规划：一是测量点规划，二是测量设备站位规划。前者是后者乃至测量实施的基础，后者是前者优化的依据。测量点规划是将待测特征离散为一组点作为理论测量点（简称测量点），用于站位规划、测量实施、模型重构、质量分析等环节，但当前测量点规划依赖于工艺人员经验，缺乏理论依据与评价机制，极易出现测量效率低、实测数据分析量大等问题。因此，合理规划测量点对测量效率、测量精度的提升具有重要意义。

Cho等［1］根据待测曲面面积、测量精度等因素利用模糊神经网络确定布点数量，为三坐标测量机布点提供了依据，但在布设过程中未考虑到自由曲面的复杂性。Lee等［2］结合Hammersley序列提出了一种布点方法，有效地减少了布点数量。宋占杰等［3］将自由曲面的曲率函数作为生成质心Voronoi结构中的密度函数，以成本函数收敛性作为算法结束准则，提出了一种全新的基于质心Voronoi的检测点采样方法。刘红军等［4］提出了实时重构布点策略，但需要在测量时对测量数据实时重构，测量效率低下。此外，诸多学者对测量点布设展开了研究，但多以特征的曲率特性作为测量点布设依据，并未考虑实际的测量不确定度［5-6］。激光雷达等测量设备在测量距离较远时，测量不确定度会达到甚至超过0.1mm，此时的测量数据难以精准反映零件的实际状态，需要采用相关算法减小测量不确定度。

本文采用非均匀有理B样条（Non-Uniform Rational B-Spline，NURBS）理论精确拟合自由曲线或曲面，并利用粒子群优化算法优化控制点和权因子提高拟合精度。对待测特征进行区域划分，包括曲率极值区和测量不确定度较大区，并针对不同的区域设计相应的布点策略确保测量点能够精确地表达待测特征。在有效保证测量效率的前提下，使测量点逆向重构特征与理论特征偏差优于10-3mm。

1 问题描述

对于复杂飞机零部件而言，曲面特征多以自由造型设计，测量点规划主要取决于工艺人员的经验。过多的测量点会延长测量周期，过少的测量点无法精确测得零件的实际状态。此外，当前测量点布设大多以曲率特性作为依据，忽视了测量不确定度对测量结果的影响，测量不确定度较大位置测量数据有时难以反映零件的实际状态。即使可以增加设备站位减小测量不确定度，但随之产生的是测量效率降低、数据转换误差累积等问题。

试验件曲率极值区域及测量不确定度较大区域分布如图1所示。图中曲率极值区域并未完全覆盖测量不确定度较大区域，在尽量避免转站前提下，仅以曲率作为布点依据难以精准测得零件的实际状态。

图1 试验件测量不确定度及曲率分布Fig.1 Measurement uncertainty and curvature distribution of test piece

2 待测特征确定性表达构建

飞机大尺寸零部件待测特征主要为曲线、曲面等基本特征，曲面特征一般是通过离散为一组曲线、采用NURBS精确拟合各条曲线进行分析。鉴于3次NURBS曲线具有C2连续性且能够满足工程中的造型需求［7］，本文的研究基于3次NURBS曲线展开。

2.1 曲线方程参数化

一条p次NURBS曲线可定义为［8］

式中：Ni，p为Bernstein基函数；u为自由变量。

自由曲线表征为NURBS初始需求解控制点Vi、权因子wi及节点矢量U，其中节点矢量通过弦长参数化确定。权因子和控制点的耦合求解过程复杂，极易出现奇异解或无解的情况，为简化计算，令wi初始值为1，控制点可由式（2）解得：

式中：Pi为型值点。

控制点求解后曲线初始方程r0（u）为

式（3）确定的拟合曲线往往不能精确表达理论曲线，一般采用增加型值点的方式提高拟合曲线的精度，但型值点数量增加到一定程度后对曲线精度的提高效果不明显，有时甚至会产生扭曲和畸变［9］，需要对wi和Vi进行优化提高曲线的拟合精度。

2.2 参数化方程优化

2.2.1 局部优化

NURBS曲线形状可通过调节控制点和权因子进行修改，权因子wi负责曲线形状的局部修改。针对拟合曲线中局部区域偏离理论曲线的情况，可通过调节wi提高曲线的拟合精度，如图2所示。

设调整后权因子为w′i，wi与w′i的关系可表示为［10］式中：d为拟合曲线与理论线间的最大偏差。

图2 权因子对曲线形状的影响Fig.2 Influence of weight factor on curve shape

对于仅有少数几个部位误差较大的情况，依据式（4）对曲线形状微调，重复采用此方法快速提高参数方程正确性。但是对于整体偏离理论曲线的情况，修改的权因子数量多，优化效率低。

2.2.2 整体优化

曲线的整体优化方法为对Vi及wi共同调整，但由于当前控制点和权因子间未形成确定的关系式，调整过程中Vi和wi的设置直接影响优化结果及速度。鉴于粒子群优化算法［11］具有收敛速度快、参数少等优点，本文采用粒子群优化算法对Vi及wi进行优化。

以一组控制点及其权因子作为一个优化粒子Zj，优化问题解的维度D＝4n（n为控制点个数，x，y，z为控制点坐标），第j个粒子定义为

设第j个粒子的最优位置为Pbestj，整个群体经历的最优位置为Gbest，第j个粒子的速度及位置更新公式分别为

式中：w为惯性权重；c1、c2为学习因子；rand1、rand2为0～1之间的随机数。设第j个粒子Zj构成的NURBS曲线为rkj（u），将rkj（u）进行离散，计算每个离散点到理论线的最小距离，形成距离集Dj＝｛dj1，dj2，…，djm｝（m为离散点个数），构造适应度函数为

式中：F为粒子个数。结合局部优化和整体优化方法，待测特征精确参数化总体流程如下：

步骤1 采用等弦长法将理论线r离散为一组型值点P＝｛P0，P1，…，Pn｝，以P为初始条件计算初始曲线方程r0（u）并生成拟合曲线。

步骤2 计算拟合曲线与理论线的最大偏差δ，若δ大于给定阈值ε0，回到步骤1增加型值点数量以减小δ，否则执行步骤3。

步骤3 判断拟合曲线与理论线偏离状态，若两曲线偏离距离δ小于阈值ε1，执行步骤7；若局部偏离执行步骤4；若整体偏离执行步骤5。

步骤4 计算r与拟合曲线最大偏差δ，按式（4）修改权因子形成新的曲线方程rk（u）。计算r与rk（u）的最大偏差δ，执行步骤5。

步骤5 若δ大于阈值ε1，重复步骤4修改下一处权因子，直到δ变化不明显或小于ε1，执行步骤7。

步骤6 采用粒子群优化算法对控制点及权因子整体优化，若δ变化不明显或小于ε1，执行步骤7。

步骤7 输出最终的控制点、权因子及参数化曲线方程r（u）。

3 待测特征测量点布设

飞机大型结构件、薄壁等零件易产生装配变形，曲率较大的位置变形较难控制［12］，且在测距较远时测量不确定度较大，测量结果不准确。针对上述情况，综合考虑曲率特性和测量不确定度等因素布设测量点。

3.1 曲线测量点布设

3.1.1 曲率极值点求解

曲率极值点是描述曲线形状的关键点，为了获取曲率极值处的变形情况，需要在曲率极值点处布设测量点。设曲线的参数方程为r（u）＝｛x（u），y（u），z（u）｝，空间曲线在某一处的曲率k（u）表示为［13］

式中：r′（u）与r″（u）分别为曲线r（u）的一阶导数与二阶导数，将其代入式（9）可得

对k（u）求一阶导数为

曲率极值点为∂k／∂u＝0的解，但此方程求解过于复杂，针对这种情况本文提出了一种基于曲率导数的曲率极值点快速定位方法，具体方法流程如下：

步骤1 将曲线方程r（u）离散，记录每个离散点对应的ui。

步骤2 计算每个离散点的曲率导数J（ui），J（ui）＝∂k／∂u｜u＝ui。

步骤3 计算相邻两点曲率导数积Di＝J（ui）·J（ui+1），执行步骤4。

步骤4 遍历所有的Di值，对Di≤0的两点T0、Q0和ui、ui+1保存，对保存的点对逐一执行步骤5和步骤6。

步骤5 取ui、ui+1的中间值um代入r（u），在T0、Q0中间生成点M0。分别计算J（um）·J（ui）与J（um）·J（ui+1）的值。若J（um）·J（ui）≤0，以中间点取代Q0，反之，以中间点取代T0。执行步骤6。

步骤6 重复步骤5更新T0、Q0，直到T0、Q0距离小于给定阈值ε2，以两点中间点作为曲率极值点。

图3为曲线某一个曲率极值点求解，左边界由T0更新至T4，右边界点由Q0更新至Q4。

图4为某零件边界曲线的曲率极值求解，圆“·”为采用本文方法求解出的曲率极值点，“”为CATIA（Computer Aided Three-dimensional Interactive Application）中曲率变化曲线上的极值点。从图中可以看出，采用本文方法求解的曲率极值与CATIA中的曲率极值基本一致。

图3 曲率极值点搜索Fig.3 Search of curvature extreme points

图4 曲线曲率极值点求解结果Fig.4 Solving result of curve curvature extreme points

3.1.2 测量不确定度评估

为降低测量不确定度对测量结果的影响，在测量不确定度较大区域加密布设测量点，通过间接平差的方式减小此影响。本文以球坐标测量系统为研究对象对测量不确定度评估。任一测量点M（x，y，z）在测量坐标系的坐标可由距离l、方位角α及天顶角β解得。

由于设备、环境因素，l、α、β采集的过程中均存在着一定的误差，为评价测量点在空间的不确定度，建立不确定度椭球模型［14-15］如图5所示。

不确定椭球的半轴长Ux′、Uy′、Uz′分别为［16-17］

式中：σα、σβ与测量系统的角度编码器分辨率有关，为仪器制造商提供的固定值；σl由最小测距不确定分量和测距不确定度距离放大系数组成。

图6为曲线轮廓度及测量不确定度允许偏差示意图，测量不确定度允许偏差范围一般为轮廓度的十分之一。图中两端及曲率较大处测量不确定度椭球较大，有的甚至超出不确定度允许偏差范围。为了减小测量不确定度，提出一种以测量不确定度为依据的测量点布设方法。具体流程如下：

步骤1 选取合适位置作为测量设备站位E（x，y，z）。

步骤2 根据曲线方程r（u），对曲线进行离散，生成离散点集G＝｛G0，G1，…，Gn｝。

图5 不确定度椭球模型Fig.5 Uncertainty ellipsoid model

图6 曲线轮廓度公差带及测量不确定度允差Fig.6 Tolerance zone and measurement uncertainty tolerance of curve

步骤3 计算每个离散点Gi到E的距离li，根据式（13）计算出每个椭球的半轴长。

步骤4 以Gi为原点，Gi与E连线为u轴，在x′-y′-z′坐标系构建单侧不确定度椭球并生成椭球离散点。

步骤5 通过坐标转换将不确定度椭球离散点由x′-y′-z′坐标系转换到x-y-z坐标系中。

步骤6 计算每个椭球离散点距离理论线的最远点，将最远点拟合成为一条曲线，对其进行光顺后作为不确定度曲线。

步骤7 在不确定度曲线超出不确定度允许范围的区域布设若干测量点。

3.2 曲面测量点布设

曲面测量点布设是将曲面离散成为一组截交线［18］，按照3.1节方法对每条截交线布设测量点。大尺寸曲面测量过程中，测量数据同样存在测量不确定度，图7为曲面轮廓度公差带和测量不确定度允差的分布情况。为了减小测量不确定度的影响，构建每个测量点的单侧不确定度椭球并生成椭球离散点，找出离散点中距离理论面最远点。将最远点拟合成为不确定度曲面，计算出不确定度曲面超出不确定度允许偏差范围的区域，并在这些区域加密布设测量点。最后将截交线离散点及根据测量不确定度增加的测量点进行筛选，对于距离过小甚至重合的两点取中间点代替。

采用上述布点策略布设的测量点有时无法精确描述理论特征，需要对测量点补充。测量点补充方法是将依据曲率及测量不确定度布设的测量点进行重构，找出重构特征与理论特征偏差较大位置，并在该位置补充测量点。以测量点重构曲面与理论曲面的最大距离偏差作为测量点能否精确描述待测特征的指标，待测特征测量点布设整体流程如图8所示，图中ε3为重构特征与理论特征的阈值。

图7 曲面轮廓度公差带及测量不确定度允差Fig.7 Tolerance zone and measurement uncertainty tolerance of surface

图8 测量点布设流程Fig.8 Process ofmeasurement points distribution

4 实验验证

为实现曲线曲面的自动布点，在CATIA［19］环境下利用CAA［20］（Components Application Architecture）技术开发测量点布设模块。软件实现包含特征预处理、特征分析及测量点规划、测量点输出3部分，如图9所示。

以图1中的试验件作为验证对象，首先提取如图10（a）所示的待测特征；随后进行特征分析和测量点布设，将曲面离散为一组截交线，对截交线按照3.1节方法规划测量点，求解出测量不确定度较大区域，并在该区域采用3.2节布点方法加密布设测量点，测量点布设结果如图10（b）所示；最后输出测量点，如图10（c）所示。为使参数方程和测量点精确拟合待测特征，阈值ε1、ε1优于10-3mm。为达到测量不确定度较大区域加密布点的目的，此区域测量点重构阈值取普通区域阈值ε3的一半。