现代飞机高机动、高气动的产品需求对装配质量优化提出了巨大挑战, 数字化测量作为飞机装配中的重要环节, 测量质量的优劣性严重影响飞机的实际装配结果。目前高精度、高效率的数字化测量技术正逐步代替传统的检测方式, 但测量实施前需基于模型进行测量规划:一是测量点规划, 二是测量设备站位规划。前者是后者乃至测量实施的基础, 后者是前者优化的依据。测量点规划是将待测特征离散为一组点作为理论测量点(简称测量点), 用于站位规划、测量实施、模型重构、质量分析等环节, 但当前测量点规划依赖于工艺人员经验, 缺乏理论依据与评价机制, 极易出现测量效率低、实测数据分析量大等问题。因此, 合理规划测量点对测量效率、测量精度的提升具有重要意义。
Cho等^[1]根据待测曲面面积、测量精度等因素利用模糊神经网络确定布点数量, 为三坐标测量机布点提供了依据, 但在布设过程中未考虑到自由曲面的复杂性。Lee等^[2]结合Hammersley序列提出了一种布点方法, 有效地减少了布点数量。宋占杰等^[3]将自由曲面的曲率函数作为生成质心Voronoi结构中的密度函数, 以成本函数收敛性作为算法结束准则, 提出了一种全新的基于质心Voronoi的检测点采样方法。刘红军等^[4]提出了实时重构布点策略, 但需要在测量时对测量数据实时重构, 测量效率低下。此外, 诸多****对测量点布设展开了研究, 但多以特征的曲率特性作为测量点布设依据, 并未考虑实际的测量不确定度^[5-6]。激光雷达等测量设备在测量距离较远时, 测量不确定度会达到甚至超过0.1 mm, 此时的测量数据难以精准反映零件的实际状态, 需要采用相关算法减小测量不确定度。
本文采用非均匀有理B样条(Non-Uniform Rational B-Spline, NURBS)理论精确拟合自由曲线或曲面, 并利用粒子群优化算法优化控制点和权因子提高拟合精度。对待测特征进行区域划分, 包括曲率极值区和测量不确定度较大区, 并针对不同的区域设计相应的布点策略确保测量点能够精确地表达待测特征。在有效保证测量效率的前提下, 使测量点逆向重构特征与理论特征偏差优于10^-3 mm。
1 问题描述 对于复杂飞机零部件而言, 曲面特征多以自由造型设计, 测量点规划主要取决于工艺人员的经验。过多的测量点会延长测量周期, 过少的测量点无法精确测得零件的实际状态。此外, 当前测量点布设大多以曲率特性作为依据, 忽视了测量不确定度对测量结果的影响, 测量不确定度较大位置测量数据有时难以反映零件的实际状态。即使可以增加设备站位减小测量不确定度, 但随之产生的是测量效率降低、数据转换误差累积等问题。
试验件曲率极值区域及测量不确定度较大区域分布如图 1所示。图中曲率极值区域并未完全覆盖测量不确定度较大区域, 在尽量避免转站前提下, 仅以曲率作为布点依据难以精准测得零件的实际状态。

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

图选项

2 待测特征确定性表达构建 飞机大尺寸零部件待测特征主要为曲线、曲面等基本特征, 曲面特征一般是通过离散为一组曲线、采用NURBS精确拟合各条曲线进行分析。鉴于3次NURBS曲线具有C²连续性且能够满足工程中的造型需求^[7], 本文的研究基于3次NURBS曲线展开。
2.1 曲线方程参数化 一条p次NURBS曲线可定义为^[8]

(1)

式中:N_{i, p}为Bernstein基函数; u为自由变量。
自由曲线表征为NURBS初始需求解控制点V_i、权因子w_i及节点矢量U, 其中节点矢量通过弦长参数化确定。权因子和控制点的耦合求解过程复杂, 极易出现奇异解或无解的情况, 为简化计算, 令w_i初始值为1, 控制点可由式(2)解得:

(2)

式中:P_i为型值点。
控制点求解后曲线初始方程r⁰(u)为

(3)

式(3)确定的拟合曲线往往不能精确表达理论曲线, 一般采用增加型值点的方式提高拟合曲线的精度, 但型值点数量增加到一定程度后对曲线精度的提高效果不明显, 有时甚至会产生扭曲和畸变^[9], 需要对w_i和V_i进行优化提高曲线的拟合精度。
2.2 参数化方程优化
2.2.1 局部优化 NURBS曲线形状可通过调节控制点和权因子进行修改, 权因子w_i负责曲线形状的局部修改。针对拟合曲线中局部区域偏离理论曲线的情况, 可通过调节w_i提高曲线的拟合精度, 如图 2所示。

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

图选项

设调整后权因子为w_i^′, w_i与w_i^′的关系可表示为^[10]

(4)

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

2.2.2 整体优化 曲线的整体优化方法为对V_i及w_i共同调整, 但由于当前控制点和权因子间未形成确定的关系式, 调整过程中V_i和w_i的设置直接影响优化结果及速度。鉴于粒子群优化算法^[11]具有收敛速度快、参数少等优点, 本文采用粒子群优化算法对V_i及w_i进行优化。
以一组控制点及其权因子作为一个优化粒子Z_j, 优化问题解的维度D=4n(n为控制点个数, x, y, z为控制点坐标), 第j个粒子定义为

(5)

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

(6)

(7)

式中:w为惯性权重; c₁、c₂为学习因子; rand₁、rand₂为0~1之间的随机数。设第j个粒子Z_j构成的NURBS曲线为r_j^k(u), 将r_j^k(u)进行离散, 计算每个离散点到理论线的最小距离, 形成距离集D_j={d_j1, d_j2, …, d_jm}(m为离散点个数), 构造适应度函数为

(8)

式中:F为粒子个数。结合局部优化和整体优化方法, 待测特征精确参数化总体流程如下:
步骤1??采用等弦长法将理论线r离散为一组型值点P={P₀, P₁, …, P_n}, 以P为初始条件计算初始曲线方程r⁰(u)并生成拟合曲线。
步骤2??计算拟合曲线与理论线的最大偏差δ, 若δ大于给定阈值ε₀, 回到步骤1增加型值点数量以减小δ, 否则执行步骤3。
步骤3??判断拟合曲线与理论线偏离状态, 若两曲线偏离距离δ小于阈值ε₁, 执行步骤7;若局部偏离执行步骤4;若整体偏离执行步骤5。
步骤4??计算r与拟合曲线最大偏差δ, 按式(4)修改权因子形成新的曲线方程r^k(u)。计算r与r^k(u)的最大偏差δ, 执行步骤5。
步骤5??若δ大于阈值ε₁, 重复步骤4修改下一处权因子, 直到δ变化不明显或小于ε₁, 执行步骤7。
步骤6??采用粒子群优化算法对控制点及权因子整体优化, 若δ变化不明显或小于ε₁, 执行步骤7。
步骤7??输出最终的控制点、权因子及参数化曲线方程r(u)。
3 待测特征测量点布设 飞机大型结构件、薄壁等零件易产生装配变形, 曲率较大的位置变形较难控制^[12], 且在测距较远时测量不确定度较大, 测量结果不准确。针对上述情况, 综合考虑曲率特性和测量不确定度等因素布设测量点。
3.1 曲线测量点布设
3.1.1 曲率极值点求解 曲率极值点是描述曲线形状的关键点, 为了获取曲率极值处的变形情况, 需要在曲率极值点处布设测量点。设曲线的参数方程为r(u)={x(u), y(u), z(u)}, 空间曲线在某一处的曲率k(u)表示为^[13]

(9)

式中:r′(u)与r″(u)分别为曲线r(u)的一阶导数与二阶导数, 将其代入式(9)可得

(10)

其中:, B=y′(u)z″(u)－y″(u)z′(u), C=x′(u)z″(u)－x″(u)z′(u), E=x′(u)y″(u)－x″(u)y′(u)。
对k(u)求一阶导数为

(11)

曲率极值点为的解, 但此方程求解过于复杂, 针对这种情况本文提出了一种基于曲率导数的曲率极值点快速定位方法, 具体方法流程如下:
步骤1??将曲线方程r(u)离散, 记录每个离散点对应的u_i。
步骤2??计算每个离散点的曲率导数J(u_i), 。
步骤3??计算相邻两点曲率导数积D_i=J(u_i)·J(u_i+1), 执行步骤4。
步骤4??遍历所有的D_i值, 对D_i≤0的两点T₀、Q₀和u_i、u_i+1保存, 对保存的点对逐一执行步骤5和步骤6。
步骤5??取u_i、u_i+1的中间值u_m代入r(u), 在T₀、Q₀中间生成点M₀。分别计算J(u_m)·J(u_i)与J(u_m)·J(u_i+1)的值。若J(u_m)·J(u_i)≤0, 以中间点取代Q₀, 反之, 以中间点取代T₀。执行步骤6。
步骤6??重复步骤5更新T₀、Q₀, 直到T₀、Q₀距离小于给定阈值ε₂, 以两点中间点作为曲率极值点。
图 3为曲线某一个曲率极值点求解, 左边界由T₀更新至T₄, 右边界点由Q₀更新至Q₄。

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

图选项

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

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

图选项

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

(12)

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

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

图选项

不确定椭球的半轴长U_x′、U_y′、U_z′分别为^[16-17]

(13)

式中:σ_α、σ_β与测量系统的角度编码器分辨率有关, 为仪器制造商提供的固定值; σ_l由最小测距不确定分量和测距不确定度距离放大系数组成。
图 6为曲线轮廓度及测量不确定度允许偏差示意图, 测量不确定度允许偏差范围一般为轮廓度的十分之一。图中两端及曲率较大处测量不确定度椭球较大, 有的甚至超出不确定度允许偏差范围。为了减小测量不确定度, 提出一种以测量不确定度为依据的测量点布设方法。具体流程如下:

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

图选项

步骤1??选取合适位置作为测量设备站位E(x, y, z)。
步骤2??根据曲线方程r(u), 对曲线进行离散, 生成离散点集G={G₀, G₁, …, G_n}。
步骤3??计算每个离散点G_i到E的距离l_i, 根据式(13)计算出每个椭球的半轴长。
步骤4??以G_i为原点, G_i与E连线为u轴, 在x′-y′-z′坐标系构建单侧不确定度椭球并生成椭球离散点。
步骤5??通过坐标转换将不确定度椭球离散点由x′-y′-z′坐标系转换到x-y-z坐标系中。
步骤6??计算每个椭球离散点距离理论线的最远点, 将最远点拟合成为一条曲线, 对其进行光顺后作为不确定度曲线。
步骤7??在不确定度曲线超出不确定度允许范围的区域布设若干测量点。
3.2 曲面测量点布设 曲面测量点布设是将曲面离散成为一组截交线^[18], 按照3.1节方法对每条截交线布设测量点。大尺寸曲面测量过程中, 测量数据同样存在测量不确定度, 图 7为曲面轮廓度公差带和测量不确定度允差的分布情况。为了减小测量不确定度的影响, 构建每个测量点的单侧不确定度椭球并生成椭球离散点, 找出离散点中距离理论面最远点。将最远点拟合成为不确定度曲面, 计算出不确定度曲面超出不确定度允许偏差范围的区域, 并在这些区域加密布设测量点。最后将截交线离散点及根据测量不确定度增加的测量点进行筛选, 对于距离过小甚至重合的两点取中间点代替。

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

图选项

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

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

图选项

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

图 9 自动布点模块架构 Fig. 9 Architecture of automatic points distribution module

图选项

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

图 10 试验件局部曲面片布点过程 Fig. 10 Points distribution process of partial surface of test piece

图选项

表 1为粒子群优化算法所采用的参数, v_max、v_min分别为最大、最小更新速度, N_p为迭代次数。
表 1 粒子群优化算法参数 Table 1 Parameters of particle swarm optimization algorithm

参数	数值
c1、c2	2.05
vmax、vmin	0.01
w	0.9
Np	200

表选项






为验证布点合理性, 采用传统布点方法在试验件上布点密度为20×20和10×10的测量点, 如图 11(a)和图 11(b)所示, 图 11(c)为采用本文方法布设的测量点。利用激光雷达MV260作为测量工具(其角度精度6.8 μm/m, 距离精度10 μm+2.5 μm/m), 对3种情况进行现场测量。

图 11 不同布点密度布点结果对比 Fig. 11 Comparison of points distribution under different measurement points density

图选项

分别构建3种布点方法中的理论点的单侧不确定度椭球, 将不确定度椭球离散, 计算每个椭球离散点集中距离理论面的最远点, 将最远点拟合为不确定度曲面, 以不确定度曲面距离理论面的最大距离作为拟合精度, 3种布点方法的测量效率及拟合精度如表 2所示。试验件测量点分布情况如图 12所示。
表 2 不同布点密度下测量结果对比 Table 2 Comparison of measurement result under different measurement points density

布点密度	布点数量	测量时长/min	重构孔洞数量	拟合精度/mm
20×20	8 868	98	21	0.037 3
本文密度	9 986	110	5	0.019 4
10×10	19 958	220	4	0.015 6

表选项

图 12 试验件测量点分布 Fig. 12 Measurement points distribution of test piece

图选项

从表 2中可以看出, 采用本文方法布点测量效率介于布点密度为10×10和20×20之间, 重构的孔洞数量、拟合精度与布点密度为10×10的结果相当, 但明显优于布点密度为20×20的测量点重构的孔洞数量及拟合精度, 充分体现了本文布点策略的合理性及可行性。
5 结论 针对航空大尺寸零部件测量点规划的基础问题, 本文分别从待测特征曲率特性、测量不确定度角度出发对测量点进行规划, 研究了测量点布设的理论依据及布设方法。主要结论如下:
1) 将待测特征精确参数化, 并对参数方程优化实现拟合精度优于10^-3 mm, 为复杂特征的曲率及测量不确定度分析提供依据。
2) 为复杂曲线曲面测量点规划提供了理论依据, 在保证测量点精确表征待测特征的条件下控制测量点数量, 使测量点逆向特征与理论特征的最大偏差优于10^-3 mm。
3) 开发了测量点自动布设模块, 提高测量点布设效率及规范性。

参考文献

[1]	CHO M W, LEE H, YOON G S, et al. A feature-based inspection planning system for coordinate measuring machines[J]. International Journal of Advanced Manufacturing Technology, 2005, 26(9-10): 1078-1087.
[2]	LEE G, MOU J, SHEN Y. Sampling strategy design for dimensional measurement of geometric features using coordinate measuring machine[J]. International Journal of Machine Tools & Manufacture, 1997, 37(7): 917-934.
[3]	宋占杰, 张美, 何改云, 等. 基于质心Voronoi结构的自由曲面布点策略[J]. 吉林大学学报(工学版), 2013, 43(1): 34-38. SONG Z J, ZHANG M, HE G Y, et al. Sculptured surface point distribution strategy based on Voronoi tessellation[J]. Journal of Jilin University(Engineering and Technology Edition), 2013, 43(1): 34-38. (in Chinese)
[4]	刘红军, 叶文静, 纪俐. 基于实时重构的自由曲面自适应布点方法[J]. 中国机械工程, 2017, 28(17): 2090-2094. LIU H J, YE W J, JI L. Adaptive distribution of inspection points for free-form curved surfaces based on real-time reconstruction[J]. China Mechanical Engineering, 2017, 28(17): 2090-2094. (in Chinese)
[5]	张虎, 张润, 于连栋. 基于CAD模型的三坐标测量机测量点分布规划[J]. 自动化与仪表, 2019, 34(2): 1-4. ZHANG H, ZHANG R, YU L D. Measurement point distribution planning of coordinate measuring machine based on CAD model[J]. Automation & Instrumentation, 2019, 34(2): 1-4. (in Chinese)
[6]	闫如忠, 张文辉. 基于离散曲率的自由曲面自适应测量技术[J]. 制造业自动化, 2018, 40(4): 153-156. YAN R Z, ZHANG W H. An adaptive measurement technology of free-form surface based on dispersed curvature[J]. Manufacturing Automation, 2018, 40(4): 153-156. (in Chinese)
[7]	魏栋.面向复杂曲面加工的NURBS曲线逼近及插补算法研究[D].杭州: 浙江大学, 2017. WEI D.Research on NURBS curve approximation and interpolation algorithm for complex surface machining[D].Hangzhou: Zhejiang University, 2017(in Chinese). http://cdmd.cnki.com.cn/Article/CDMD-10335-1017047246.htm
[8]	朱心雄. 自由曲线曲面造型技术[M]. 北京: 科学出版社, 2000: 152-160. ZHU X X. Free curve and surface modeling technology[M]. Beijing: Science Press, 2000: 152-160. (in Chinese)
[9]	JAUCH J, BLEIMUND F, RHODE S, et al. Recursive B-spline approximation using the Kalman filter[J]. Engineering Science and Technology, an International Journal, 2017, 20(1): 28-34.
[10]	闫艳.NURBS曲线形状修改的研究[D].郑州: 郑州大学, 2012. YAN Y.The research on the modification of NURBS curves[D].Zhengzhou: Zhengzhou University, 2012(in Chinese).
[11]	KENNEDY J, EBERHART R.Particle swarm optimization[C]//Proceedings of ICNN'95-International Conference on Neural Networks.Piscataway: IEEE Press, 2002: 1942-1948.
[12]	韩志仁, 王刚, 彩辉, 等. 基于CATIA/CAA的复杂边界平面和曲面自动布点技术研究[J]. 制造业自动化, 2017, 39(9): 59-63. HAN Z R, WANG G, CAI H, et al. Research on generating point automaticly in complex boundary plane and surface based on CATIA/CAA[J]. Manufacturing Automation, 2017, 39(9): 59-63. (in Chinese)
[13]	陈维桓. 微分几何[M]. 北京: 北京大学出版社, 2006: 23-47. CHEN W H. Differential geometry[M]. Beijing: Peking University Perss, 2006: 23-47. (in Chinese)
[14]	黄鹏, 王青, 李江雄, 等. 激光跟踪仪三维坐标转换综合优化方法[J]. 计算机集成制造系统, 2015, 21(11): 2912-2920. HUANG P, WANG Q, LI J X, et al. Comprehensive optimization for three-dimension coordinate transformation of laster tracker[J]. Computer Intergrated Manufacturing System, 2015, 21(11): 2912-2920. (in Chinese)
[15]	PREDMORE C R. Bundle adjustment of multi-position measurements using the Mahalanobis distance[J]. Precision Engineering, 2010, 34(1): 113-123.
[16]	DENG Z, LI S, HUANG X. Uncertainties evaluation of coordinate transformation parameters in the large-scale measurement for aircraft assembly[J]. Sensor Review, 2018, 38(4): 542-550.
[17]	DENG Z, LI S, HUANG X. Coordinate transformation uncertainty analysis and reduction using hybrid reference system for aircraft assembly[J]. Assembly Automation, 2018, 38(4): 487-496.
[18]	ZABOTIN I, KAZAEVA K E. The procedure for the complete updating of immersive sets in one cutting plane method[J]. Journal of Physics:Conference Series, 2019, 1158: 042041.
[19]	张俐, 江春, 李承文. CATIA平台下的机身数字化对接测量软件开发与应用[J]. 制造业自动化, 2017, 39(2): 129-133. ZHANG L, JIANG C, LI C W. Measurement software development and application of fuselage docking with CATIA platform[J]. Manufacturing Automation, 2017, 39(2): 129-133. (in Chinese)
[20]	张辉, 李泷杲. 基于CATIA二次开发的制动盘摩擦半径计算技术研究[J]. 制造业自动化, 2018, 40(2): 74-76. ZHANG H, LI S G. Research on friction radius calculation of brake disc based on secondary development of CATIA[J]. Manufacturing Automation, 2018, 40(2): 74-76. (in Chinese)