机器人原位铣削技术应用越来越广泛，三维数字化测量在保证飞机装配蒙皮铣削精度方面起到了重要作用。如采用激光测距、相位光栅、三角定位以及立体视觉等原理的扫描设备，与机器人协同完成自动铣削/制孔的部件型面扫描、全机外形测量等工作，实现了点云数据的快速采集，显著提高了飞机装配精度。但对于飞机大型曲面（面积超过20 m2），一次测量点云数量数以亿计，直接处理点云数据极大增加计算机运算负担，占用大量存储容量和内存空间，而且提取用于质量评价或加工的有效信息需耗费大量时间，显著降低了产品检测和自动化铣削的效率。因此，为提高测量数据处理效率，对点云进行压缩处理，开展面向机器人原位铣削的数据快速精简技术研究势在必行。
随着数字化测量技术的发展，非接触光学测量设备，如激光跟踪仪、激光雷达、ATOS、Creaform 扫描仪等获得大量应用，测量得到的原始点云数据存在规模大、测点噪音与层叠等问题，严重影响数据分析精度和效率。为此大量学者对点云数据的滤波和精简进行研究，开发了多种处理算法[1]，国外学者率先开展了散乱点云滤波方法研究，常用方法有高斯滤波[1]、双边滤波[2–4]、拉普拉斯滤波[5]等，高斯滤波通过加权平均处理抑制噪音影响，但不适用于高曲率和薄壁特征的降噪。拉普拉斯滤波计算简单，适应平滑曲面点云光顺，对尖锐和小曲率半径区域点云有损伤。双边滤波应用较多，采用邻域点的加权平均进行处理，有效保留了几何边缘特征，适用于孔洞、筋条、空间曲线边界的处理。国内学者基于前人成果，开展了大量点云精简算法研究，主要包括包围盒法[6]、随机采样法[7]、角度法[8]、均匀网格法[9]和曲率采样法[10]，包围盒法基于八叉树理论，将所有点云数据包围在内并进行有序分裂，通过采样点的融合进行点云压缩，有着良好的序列性和分层性，便于空间划分，提高了点云轻量化速度。随机采样法易于实现，按设定精简率对点云数据进行随机索引删除，精简速度快，但容易造成细节特征的损失，应用较少。角度法计算采样点与邻近两点连线间的夹角，并与设定的角度阈值比较，控制点云删除与保留，完成点云精简，计算效率较高，但造成曲率变化小区域精度损失。均匀网格法将点云分配给与扫描方向垂直的均匀平面方格，计算同一方格内各点到方格平面的距离并排序，压缩掉除“中值点”以外测点，从而完成压缩，该法对尖锐和小曲率半径区域特征损伤大。曲率采样法对点云数据进行曲面拟合，得到各处曲率估算值，曲率较小区域采用大压缩率，曲率较大区域采用较小压缩率，该方法对曲面特征表述有利。基于上述理论开发的商业软件有Geomagic、Polyworks 等，基本满足了用户使用需求，但商业软件需要大量的人机交互操作，不便于自动化设备的集成化应用。为此，国内学者针对飞机质量控制和自动加工需求，开发了快速检测分析系统[11]，将数字化测量技术与自动化加工设备集成[12]，在蒙皮下陷检测和加工方面取得了较好的应用效果，保证了检测精度，提高了检测分析效率。
综上所述，为满足飞机装配机器人原位铣削快速检测分析需求，亟须开展离散点云数据快速精简技术研究，保证关键特征的检测精度，实现大型构件的高效检测与分析，提高机器人铣削自动化程度，促进飞机装配工艺技术水平的不断完善。
1 数据精简算法
飞机装配原位铣削对关键特征测量的时效性要求高，通过测量点云降噪和压缩处理，保留关键特征信息，降低点云数据规模，为铣削边界信息的快速提取奠定基础。
1.1 双边滤波
双边滤波降噪效果明显，且能够保持点模型的几何特征，在扫描图像处理中具有良好的鲁棒性，其核心思想是在点云数据上设置多个邻域，并依据判定准则为域内点赋予差异化的滤波权值，进而实现点云的滤波降噪。飞机原位铣削关键特征包括圆角、筋条、加工边界等，以下采用双边滤波法处理点云数据，可有效保留特征边缘特征，降低滤波对边界特征的损伤。
二维双边滤波根据图像领域内灰度值进行处理，但不适于点云数据滤波。本研究对二维滤波进行改进，首先根据点云数据拟合确定表面法向，再以法向为参考计算出领域内各点高度值，然后仿照灰度值处理制定三维点云双边滤波算法，具体为
其中，
式（1）和（2）中，p′代表点p 经过双边滤波处理后的计算值；p为处理点；q为p 邻域内点；S为点云域；np为单位法向量；σd 和σr 分别代表定义域和值域中高斯滤波器的滤波器核，决定了滤波过程中邻域内点云值的权重。定义域滤波是空间邻域内处理点的加权平均值，其中q 和p 之间距离越大，相应权重系数就越小，反之则越大。类似的，值域滤波也是一个加权平均值，对处理点数值邻域内的点进行处理，其中坐标差值越大，权重系数就越小，反之则越大。
Wp 是一个标准量，其值为
是空间函数，定义为
是范围函数，定义为
式（4）中，d（p，q）为两点间距离差；式（5）中，δ（p，q）为两点间沿法向高度差值。
1.2 点云精简
点云精简的方法较多，工程上常用的有均匀精简法和曲率适应性精简法，均匀精简主要用于相对平缓（曲率变化小）的曲面，曲率适应性精简主要用于曲率变化大具有高曲率特征和薄壁特征的曲面。点云精简算法大都基于八叉树理论制定，采用八叉树划分，可处理任意拓扑结构点云数据，便于数据的组织和压缩处理。八叉树划分的压缩算法主要流程分八叉树构建和数据编码两个阶段。
首先根据点云数据极限值建立包围盒，保证所有数据均处于包围盒内。已知点云数据集P = {p1，p2，…，pm}，计算点云集最大坐标值（xmax，ymax，zmax）和最小坐标值（xmin，ymin，zmin），建立边长为lx、ly、lz 的包围盒，有
然后对包围盒进行均匀划分，一个栅格被均分为8 个子栅格。图1为栅格划分和八叉树结构的示例，栅格细化的过程中，判断栅格是否为空，只有非空栅格才会进一步细化，细化程度根据精度要求进行。
图1 八叉树构建
Fig.1 Octree building
如图2所示，按照文献[1]，为了避免部分测点处于包围盒边界，适当增加包围盒的3 条边lx、ly、lz 边长，以保证所有点云数据点均位于包围盒内。根据点云规模和精简比例要求，将包围盒划分为s 个边长为a 的小栅格。
图2 点云栅格精简
Fig.2 Voxel-grid simplificatio
式中， 表示向上取整运算。令点云包含测量点数为m，定义每个栅格中含有测量点数均值为
，则有
，可以得出每个细化栅格边长估计值为
定义点云集为P，通过八叉树划分，将所有点云观测点pi∈P 投射到相应的栅格中，每个栅格其内的点被压缩，最后观测点pi 所在包围盒栅格空间位置[gix giy giz]T 表示为
由式（9）即得压缩后采样点。通过如此划分，八叉树中每一个层次上的非空栅格集合就成为点云的一个近似表示，各层的近似模型就构成了原始点云的一个多细节层次的表示。栅格的每次划分均对栅格进行位置编码，非空为1，反之为0，这样对于每次栅格细分，用一个固定的顺序遍历所有子栅格，得到占位编码信息，而且与之关联的位置和法线信息均被编码。为了降低占位编码的熵，将占位码重新排序，把非空集都推到前面，便于后续遍历索引。
兼顾点云压缩规模和模型标定精度，在曲面平坦区域按上述算法保留均匀测点，在曲面高曲率区域采用基于局部法矢量变化和平均曲率的适应性精简算法，保留相对密集测点。假设某一栅格内各观测点为pi（i=1，…，m），则栅格内对应观测点平均曲率为ρi（i=1，…，m），可得该栅格平均曲率为
根据精度要求制定曲率精简阈值δρ，若
，则仅保留栅格内曲率最接近ρ-的观测点；若
，则对该栅格再行细分，细分栅格边长a′为
式中，a为原栅格长度；θmax为原栅格内任意两观测点法矢最大夹角；δθ为给定的法矢变化阈值； 代表向下取整运算符号。由于各栅格最大法矢夹角θmax 不尽相同，细分栅格长度与平均曲率
亦不相同，若
则接着细分，直至满足
，保留栅格内曲率最接近
的采样点。
2 点云精简验证与集成应用
基于上述算法，通过科研课题合作方式，与华中科技大学联合开发应用程序iPoint3D-AeroSpect，如图3所示，采用典型样件对点云数据精简效果进行验证，并与自动铣削设备进行集成应用，证明本文方法的有效性。
图3 iPoint3D-AeroSpect 软件界面
Fig.3 Software interface of iPoint3D-AeroSpect
2.1 测量点云精简验证
飞机表面对缝处多采用锯齿结构以提高隐身性，该类零件具有典型边界特征，扫描测量与特征提取直接影响自动化铣削精度。仿照锯齿口盖类零件（长1102 mm，宽695 mm），采用高精度机床加工试验样件，模拟实际装机零件加工后喷漆处理，漆层厚度不超过10 μm，同时喷漆可降低金属零件反光性，便于点云数据采集。试验样件准备完成后，采用三坐标测量机（测量精度3 μm+4L/1000 μm，L为长度）进行精确标定，作为验证检测标准。
采用三维扫描仪进行数据采集（单幅测量误差不超过0.03 mm），能够获得精细特征细节信息，原始点云数据如图4（a）所示（其中圆孔为夹持模块，用于自动铣削夹紧和减振），点云量为119970312 个，包括大量噪点。采用iPoint3D-AeroSpect 应用程序，通过双边滤波进行点云数据降噪处理，并对点云数据进行压缩，最大限度保留边界特征信息，处理后零件边缘区域点云数据如图4（b）所示，点云数量为18892962 个，显著降低了点云规模，压缩比达到84.25%。
图4 点云数据精简
Fig.4 Simplification of point clou
对上述处理后点云数据进行详细分析，如图5所示，提取轮廓边界关键特征点，然后根据曲线特征拟合进行均匀精简，最后采用Nurbs 样条曲线进行拟合，得到试验样件边界特征的数学描述。
图5 轮廓边界提取
Fig.5 Extraction boundary of contour
装配配合间隙及均匀性对飞机隐身性影响显著，影响因素包括总体轮廓、边夹角以及圆角特征，基于得到边界特征曲线，重点对上述特征进行验证。将三坐标测量机对上述特征的精确测量作为数据参考，对比结果见表1。
按试验样件的边界特征种类划分，与三坐标测量机数据对比，表1所列数据为各特征误差最大值，其中轮廓误差0.05 mm、角度误差0.021°、圆弧角半径误差0.03 mm，满足边界质量评价和自动化应用需求。综上所述，本文研究算法实现了点云数据的高效精简，有效保留了边界信息，降低了点云压缩算法对边界特征测量精度的损失，特征边界精度满足装配蒙皮铣削要求。
表1 处理结果对比
Table 1 Contrast of processing results
最大误差项圆角半径/mm标定值拟合值偏差2.4922.5220.03边夹角（锐角）/（°）75.00175.0220.021边夹角（钝角）/（°）141.998142.0170.019总体轮廓/mm——0.05
2.2 机器人铣削集成应用
铣削工艺流程如图6所示，由于飞机零部件具有大型、薄壁弱刚性、形状复杂等特点，为消除误差累积、零件自身变形等不确定性因素影响，机器人铣削系统需要现场采集数据并导入应用程序iPoint3D-AeroSpect，进行数据精简，然后通过数据匹配与飞机坐标系对齐，按照测量与铣削工艺要求进行特征区域分割并提取特征，明确待加工产品的实际状态，以此为基础完成加工路径规划和程序生成，程序下发至机器人铣削系统完成加工，最后通过测量评价加工质量。
图6 铣削工艺流程图
Fig.6 Milling process flow chart
软件所有功能接口开放，铣削过程中控制系统可任意调用软件功能函数，实现了同机器人铣削设备高效集成，如图7所示，显著降低了人机交互操作频率，并通过离散点云数据精简处理，显著提高了加工特征信息提取的效率，为铣削加工提供及时高精度的数据支持，避免了质量偏差，满足了自动铣削设备的时效性需求。以某型部件为例，初始测量点云数据203352569 个，应用本文技术前，平均单件数据处理、分析与加工数据提取需要熟练操作人员点击鼠标不少于200 次，耗时不低于28 min。应用本文技术后，所有数据处理过程无需人工干预，点云规模压缩为36603462 个，平均单件数据处理、分析与加工数据提取耗时降为7 min，显著提高了铣削效率。
图7 机器人自动测量–铣削系统集成应用
Fig.7 Integration application of robot automatic measuring–milling system
3 结论
（1）针对飞机装配蒙皮修切与连接过程对数字化测量技术的需求，开发了适应三维点云的双边滤波处理方法，实现了空间点云数据的降噪处理，降低滤波对边界特征的损伤，有效保留了复杂构件圆角、斜边、空间曲线边界等关键特征。
（2）基于八叉树理论，结合均匀压缩和曲率适应性压缩，开展了适应飞机装配特点的点云数据快速精简算法研究，便于数据的组织和压缩处理，有效降低点云数据规模，保证了复杂型面检测精度，为后续边界信息的快速精确提取奠定了基础。
（3）开发了点云采集和处理应用软件，以典型样件为测试对象，通过与三坐标测量结果对比，验证了本文研究的有效性，实现了点云数据的高效精简，降低了点云压缩算法对边界特征损伤，有效保留了边界信息。软件接口开放性高，实现了与自动化铣削设备的集成应用，为修正加工程序提供及时高精度的数据支持，避免了质量偏差，降低了装配返工率。
[1]李文龙,李中伟,毛金城.iPoint3D曲面检测软件开发与工程应用综述[J].机械工程学报,2020,56(7):127–150.LI Wenlong,LI Zhongwei,MAO Jincheng.The development and application review of iPoint3D software for surface inspection[J].Journal of Mechanical Engineering,2020,56(7):127–150.
[2]F L E I S H M A N S,D R O R I I,COHENOR D.Bilateral mesh denoising[J].ACM Transactions on Graphics,2003,22(3):950–953.
[3]JONES T R,DURAND F,DESBRUN M.Non-iterative,feature-preserving mesh smoothing[J].ACM Transactions on Graphics,2003,22(3):943–949.
[4]SMITH S M.SUSAN—A new approach to low level image processing[J].International Journal of Computer Vision,1997,23(1):45–78.
[5]CAO J J,TAGLIASACCHI A,OLSON M,et al.Point cloud skeletons via Laplacian based contraction[C]//2010 Shape Modeling International Conference.Aix-en-Provence:IEEE,2010.
[6]陈达枭,蔡勇,张建生.散乱点云精简的一种改进算法[J].计算机应用研究,2016,33(9):2841–2843.CHEN Daxiao,CAI Yong,ZHANG Jiansheng.Improved algorithm of simplifying scattered point cloud data[J].Application Research of Computers,2016,33(9):2841–2843.
[7]TAGLIASACCHI A,ZHANG H,COHEN-OR D.Curve skeleton extraction from incomplete point cloud[J].ACM Transactions on Graphics,2009,28(3):1–9.
[8]李新,王细洋.基于扫描线的三维激光扫描数据精简算法[J].机械设计,2010,27(4):17–20.LI Xin,WANG Xiyang.Algorithm of reduction for 3D laser scanning data based on scanning line[J].Journal of Machine Design,2010,27(4):17–20.
[9]金涛,陈建良,童水光.逆向工程技术研究进展[J].中国机械工程,2002,13(16):1430–1436.JIN Tao,CHEN Jianliang,TONG Shuiguang.A review of reverse engineering technology[J].China Mechanical Engineering,2002,13(16):1430–1436.
[10]赵龙,韦群.序列图像三维重构中点云精简算法的研究与改进[J].计算机工程与应用,2016,52(8):211–216,228.ZHAO Long,WEI Qun.Research and improvement of point clouds reducing algorithm based on a set of images in progress of 3D-reconstruction[J].Computer Engineering and Applications,2016,52(8):211–216,228.
[11]李栋,李泷杲,李琦,等.基于激光扫描的飞机蒙皮下陷特征点提取[J].航空制造技术,2018,61(9):96–101.LI Dong,LI Shuanggao,LI Qi,et al.Feature points extraction of aircraft skin sag based on laser scaning[J].Aeronautical Manufacturing Technology,2018,61(9):96–101.
[12]康晓军,向兵飞,李迎光,等.基于快速扫描的飞机大型蒙皮自适应加工技术[J].航空制造技术,2019,62(6):85–90,94.KANG Xiaojun,XIANG Bingfei,LI Yingguang,et al.Adaptive machining technology of large aircraft skin based on rapid scanning[J].Aeronautical Manufacturing Technology,2019,62(6):85–90,94.
