涡轮叶片是航空发动机关键的热端部件，其加工质量决定了整机的性能与寿命。涡轮叶片制造工艺复杂，其精铸精度要求高（如某型号叶片型面轮廓度要求0.2mm，[–0.1 mm，+0.1 mm]），首先通过镍基高温合金等材料精铸而成，随后在快速、精确配准定位的前提下，进行数控加工、特种加工等工序。其中，配准定位是通过寻找涡轮叶片测量点云与CAD 模型间最佳相对位姿实现的。
配准定位方面应用最为广泛的是ICP 算法及其改进之后的算法[1–2]，其关键步骤中，预对齐和匹配点集搜索关系到ICP 算法求解的精度和效率。预对齐算法可为ICP配准提供较好的初始位姿，避免陷入局部最优。目前，常见的预对齐方法有重心重合法[3]、6 点对齐法[4]、标签法[5]主成分分析法[6]、特征提取法[5]等。而匹配点集搜索则是为了提高ICP 配准效率，Li 等[7]使用三角网格搜索匹配点集，通过减少搜索最近点集的时间来提高算法效率。Geng 等[8]则采用了基于邻域约束的空间投影法。程军等[9]构建kd–tree 将点云空间划分为多个子空间，并利用并行计算加速查询匹配点集。
为此，本研究针对涡轮叶片配准定位的实际需求，分析了预对齐和匹配点集搜索的重要性；在自主研发CAD/CAM 软件平台上以Visual studio 2015 为开发工具，使用基于主成分分析、手工关联两种预对齐算法实现测量点云与CAD 模型具有较好的初始位姿；对网格法、赋范空间投影法、基于kd–tree 快速搜索法3 种匹配点集搜索算法进行效率对比，并确定了最优的涡轮叶片配准定位方案。
1 ICP 配准定位算法分析
经典的ICP 算法[10]是由McKay 和Besl 提出的，其核心思想是迭代调姿使距离偏差最小化，目标函数如式（1）所示，利用奇异值分解法[11]和四元组法[12]找出测量点云与CAD 模型匹配点集的空间变换矩阵，多次迭代直至目标函数满足一定的精度为止。
式中，pi（i=0，1，…，n）为待配准的测量点；qi（i=0，1，…，n）为pi 在CAD 模型上的匹配点；（qi，pi）为匹配点集；R、T 分别为待求的旋转矩阵与平移矩阵。
针对涡轮叶片实际应用情况，通常测量点云依靠光学扫描设备获取，且初始点云位姿与涡轮叶片CAD 模型相差悬殊，直接对两者进行配准定位，往往会陷入局部最优，并且效率较低。为此，在使用ICP 算法时提出如下需求。
（1）避免初始位置依赖性。
初始位置对ICP 算法的求解结果影响很大，其收敛速度及收敛解直接受制于求取的匹配点集是否精确。为此，在第2 节采用基于主成分分析预对齐与手工关联预对齐方式保证测量点云与叶片CAD 模型间较好的初始位姿。
（2）避免匹配点计算效率低。
ICP 算法迭代过程中需要频繁计算测量点与CAD模型间的最近点，即匹配点集搜索，搜索效率决定了配准定位的计算效率。为此，在第3 节对3 种点到曲面最近点的算法进行比较，包括网格法、赋范空间投影法和基于kd–tree 快速搜索法。
如图1 所示，以典型涡轮叶片光学测量点云为数据源，以叶片CAD 模型为目标源进行上述配准定位中预对齐和匹配点集搜索算法的应用研究。
图1 涡轮叶片测量点云与CAD 模型
Fig.1 Turbine blade measurement point cloud and CAD model
2 预对齐方式对比分析
开展ICP 精确配准定位最基础的是预对齐，为了便于后续匹配点集的快速搜索，测量点云尽量逼近CAD模型。
2.1 手工关联预对齐
手工预对齐方法是通过手工建立多组点–对应面的关联关系，具体流程如下：
（1） 通过交互界面，手工建立测量点云P ={p1，…，pi，…，pm }和CAD 模型面片F ={f1，…，fj，…，ft}的关联关系Mij =（pi，fj）；
（2）计算各个关联Mij 中（pi，fj）形成的匹配点Q ={q1，…，qi，…，qm}；
（3）使用ICP 算法配准对齐（P，Q）。
执行过程如图2 所示，相应结果在2.3 节中进行对比分析。
图2 手工预对齐
Fig.2 Manual pre-alignment
2.2 基于主成分分析的预对齐
基于主成分分析 （PCA）的预对齐是通过计算测量点云的协方差矩阵，求得该协方差矩阵特征值所对应的特征向量，以两两垂直的3 个特征向量为空间坐标系的3 个坐标轴，以点云数据的重心为坐标原点，建立该点云的特征坐标系。通过将测量点云和CAD 模型离散点云的特征坐标系调整到一致来达到快速预对齐的目的。
具体实施步骤：
（1）对CAD 模型进行离散，得到离散点云Q ={q1，…，qj，…，qn}。
（2）计算P 和Q 的重心坐标gp 与gq，分别对P 和Q 内各点进行去均值标准化处理。
pi = pi – gp（i = 1，2，…，m）
qj = qj – gq（j = 1，2，…，n）
（3）点云依坐标值进行矩阵转置为PT = [p1，…，pm]T =[Xp，Yp，Zp ]，QT = [q1，…，qn]T = [Xq，Yq，Zq]，分别求PT、QT 的协方差矩阵，即
（4） 计算covp、covq 特征向量[α1，α2，α3 ]和[β1，β2，β3]，并以gp 与gq 为原点建立特征坐标系Wp 与Wq。
（5）经过旋转矩阵R 和平移向量T 使得Wp 与Wq重合，从而实现基于主成分分析的预对齐。
执行过程如图3 所示，相应结果在2.3 节中进行对比分析。
图3 基于主成分分析的预对齐
Fig.3 Pre-alignment based on principal component analysis
2.3 预对齐方式对比分析
在2.1 与2.2 节中针对初始位姿差异较大的测量点云与CAD 模型，详细介绍了手工关联和基于主成分分析的预对齐方法，相应配准定位精度结果如表1 所示。可知，基于本次手工选点所得到的配准手工关联对齐结果较主成分分析对齐结果理想，但人工参与后降低了预对齐的效率。手工关联预对齐的各个关联形成匹配点的计算并不复杂，可以较快进行对齐配准，而主成分分析是计算测量点云的协方差矩阵，计算过程复杂，因此手工关联预对齐相较于主成分法耗时短。
表1 预对齐各轴偏移结果较理论值差异对比
Table 1 Comparison of deviation results of pre-aligned axes with theoretical values
基于主成分分析的预对齐绕X 轴旋转/（°）0.1290.1630.083绕Y 轴旋转/（°）0.2090.1870.154绕Z 轴旋转/（°）–45.045–46.532–42.573沿X 轴平移/mm3.2272.8402.876沿Y 轴平移/mm–6.783–6.125–5.948沿Z 轴平移/mm–0.137–0.143–0.0632对齐结果最终变换手工关联预对齐
3 匹配点集搜索算法对比分析
因叶片后续工艺需求，经过预对齐之后，采用ICP算法进一步精确配准定位叶片自由曲面CAD 模型及其相应测量点云的空间相对位姿。其中，匹配点集的计算效率问题是ICP 算法的主要问题。
3.1 网格法搜索匹配点集
网格法求取测量点到自由曲面最近点的基本思路是将自由曲面的参数域D ={（u，v）| umin≤u≤umax，vmin≤v≤vmax}分别沿着u 向和v 向等分成m 份和n 份，参数域上的m×n 个网格中心点参数值为Eij = {（ui，vj），i = 1，2，…，m；j = 1，2，…，n}，自由曲面上所对应的三维坐标点为qij，选取其中最小的距离作为最短距离，并对该中心点所对应的参数域网格继续细化，重复以上网格划分求取最近点的过程，直到前后两次最近点的差值小于一定的迭代精度终止。网格法细分曲面，求取最近点的原理如图4 所示。
图4 网格法计算匹配点原理图
Fig.4 Schematic diagram of calculating matching points by grid method
3.2 赋范空间投影法搜索匹配点集
赋范空间投影法是在充分考虑曲面的局部微分性质的基础上，利用迭代逼近的方法求取空间点到自由曲面的极值距离。但是，在复杂曲面边界处或者封闭曲面上进行最近点搜索匹配时，该方法极易不收敛。为此，对该方法进行改进：在不考虑自由曲面边界的前提下，使用迭代法求解空间点在曲面上的垂足点，判定垂足点是否在曲面的位置。若在曲面内，那么垂足点就是最近点；若不在曲面内，则根据垂足点Q 的uv 参数进行处理分析，使得最近点落在相应的曲面边界上，其基本原理如图5 所示，将最终的迭代最近点Q（u，v）后置到了Q'（u'，v'）处。
图5 赋范空间投影法边界处理示意图
Fig.5 Boundary treatment diagram of normed space projection method
在图5 中，F 为自由曲面，P 为空间点，Q（u0，v0）为F 上的迭代初始点，向量
，空间向量n、dv、du分别是Q（u0，v0）处的法矢量、u 向切矢和v 向切矢，在Q（u0，v0）处形成了放射坐标系框架，Q（u，v）为P 在放射坐标系dv、du 面上的投影点。设 （Δu，Δv）为参数（u，v）与参数 （u0，v0）的差值，l 为向量r 在法矢量n 上的投影，则
式 （2）两侧分别乘以du、dv，且由于n ·du= 0，n ·dv = 0，令r·du=A，r·dv= B，du·du= E，du·dv = F，dv·dv = G，则式 （2）变换成
du、dv 是曲面F 在初始迭代点Q（u0，v0）处的局部微分信息， 
，至此，可以顺利求出未知量Δu 和Δv，从而确定投影点Q（u+Δu，v+Δv）的坐标，再以Q（u+Δu，v+Δv）作为初始迭代点，重复执行以上过程，直到 （Δu）2+（Δv）2≤ε 时，迭代终止，垂足点Q 就是空间点P 到曲面F 的最近点。而曲面都是有边界的，所以垂足点Q 有可能不落在曲面F 内部。当垂足点Q 落在F内部时，Q 就是所求的最近点，| PQ |就是最短距离；当垂足点Q 不在F 内部时，最近点不再是垂足点Q，而应该在自由曲面边界上。
3.3 基于kd–tree 快速搜索匹配点集
kd–tree 搜索算法是一种高维空间中最近邻与近似最近邻的快速查找技术，它将整个高维空间按照一定的原则、递归地进行划分，然后查询点根据kd–tree 建立的规则，在特定空间内进行查询与回溯操作，从而找到距离函数最小的匹配点。基于kd–tree 快速搜索匹配点集算法分为两大部分，即kd–tree 的建立和基于kd–tree 最近点的快速搜索。
创建kd–tree 流程如图6 所示，具体步骤如下。
图6 创建kd–tree流程图
Fig.6 Flowchart of creating kd–tree
（1）在K 维数据集合M 中计算并选出最大方差的维度SplitDim，而后对该维度SplitDim 上的数据进行升序排列，取其中值MidVal 对该数据集合M 进行划分，得到两个子集合Mleft 和Mright；同时创建一个树节点KDTreeNode 来存储中值MidVal 所对应的高维数据点。
（2）递归地对两个子集合重复步骤（1）的过程，直至所有子集合都不能再划分为止。
而快速搜索的基本思路是基于kd–tree 这种平衡排序树的数据结构，在kd–tree 上进行比较分析确定搜索路径，从而锁定查询的极小空间区域。随后进行回溯操作以确定最近点。因此kd–tree 搜索算法计算效率较好。算法的具体实施步骤如下：
（1）在kd–tree 上进行二分查找，形成搜索路径；
（2）用搜索路径中的最后一个节点（叶节点）作为当前最近点；
（3）回溯查找，沿着搜索路径由下而上进行回溯搜索，更新最近点；
（4）当查询点的邻域与当前节点的分割超平面相交时，需要用当前节点另一侧的子分支，跳转到步骤（1），补充新的搜索路径。
3.4 匹配点计算效率的比较
对图1 所示算例进行上述3 种匹配点集的搜索，情况如表2 所示，可知，网格法的计算效率最低，同等条件下，计算耗时约为赋范空间投影法的3 倍；基于kd–tree搜索算法的效率最高，并且测试点数越多效果越明显。
表2 匹配点集计算时间对比
Table 2 Comparison of calculation time of matching point set
点数方法耗时/s网格法赋范空间投影法基于kd–tree 搜索法1470.210 0.078 0.078 5390.880 0.281 0.109 45115.100 2.294 0.608
由此，采用手工关联预对齐、基于kd–tree 快速搜索匹配点集和ICP 方法对图1 所示算例进行配准定位，在自主研发CAD/CAM 软件平台上定位效果如图7 所示。实例测量点云与CAD 模型差异区间为[–0.06 mm，+0.1 mm]，符合涡轮叶片精铸精度情况。
图7 涡轮叶片定位前与最优定位效果
Fig.7 Before turbine blade positioning and optimal positioning effect
4 结论
本文研究了涡轮叶片精铸后测量点云与CAD 模型的配准定位技术，重点对比分析了影响配准精度的预对齐方法，以及影响配准效率的匹配点集搜索方法，得到如下结论。
（1）分析了利用ICP 算法进行配准定位的特点和难点问题，指出了提高配准精度的预对齐以及提高配准效率的匹配点集搜索需求。
（2）针对初始位姿不好的测量点云，实现了基于主成分分析的预对齐和手工关联预对齐，对比分析表明手工关联预对齐具有更好的配准精度，但人工参与不可避免带来预对齐效率的降低。
（3）针对匹配点计算效率较低的问题，实现3 种匹配点的计算方法，通过比较分析，基于kd–tree 搜索算法计算效率较好。在此基础上，确定了手工关联预对齐、基于kd–tree 快速搜索匹配点集和ICP 方法的配准方案。
（4）优化出了一种配准算法组合。配准完成之后使得整体型面的余量满足后续加工要求，配准定位后可以保证用户指定的余量要求。
[1]朱琛琛. 基于ICP 算法的点云配准研究[D]. 郑州: 郑州大学, 2019.
ZHU Chenchen. Point cloud registration based on ICP algorithm research[D]. Zhengzhou: Zhengzhou University, 2019.
[2]LI X D, LI W, JIANG H Z, et al. Automatic evaluation of machining allowance of precision castings based on plane features from 3D point cloud[J]. Computers in Industry, 2013, 64(9): 1129–1137.
[3]李运川, 王晓红. 一种用于点云配准的改进迭代最近点算法[J]. 软件导刊, 2020, 19(9): 175–179.
LI Yunchuan, WANG Xiaohong. An improved iterative closest point algorithm for point cloud registration[J]. Software Guide, 2020, 19(9):175–179.
[4]李聪波, 肖卫洪, 杜彦斌, 等. 基于改进ICP 算法的损伤零部件精确配准方法[J]. 计算机集成制造系统, 2016, 22(4): 1021–1028.
LI Congbo, XIAO Weihong, DU Yanbin, et al. Fine registration method for defective parts based on improved ICP algorithm[J].Computer Integrated Manufacturing Systems, 2016, 22(4): 1021–1028.
[5]严思杰, 周云飞, 彭芳瑜, 等. 大型复杂曲面加工工件定位问题研究[J]. 中国机械工程, 2003, 14(9): 737–740.
YAN Sijie, ZHOU Yunfei, PENG Fangyu, et al. Research on localization of the workpieces with large sculptured surfaces[J]. China Mechanical Engineering, 2003, 14(9): 737–740.
[6]刘哲, 周天, 彭东东, 等. 一种改进的基于PCA 的ICP 点云配准算法研究[J]. 黑龙江大学自然科学学报, 2019, 36(4): 473–478, 505.
LIU Zhe, ZHOU Tian, PENG Dongdong, et al. An improved ICP point cloud registration algorithm based on PCA[J]. Journal of Natural Science of Heilongjiang University, 2019, 36(4): 473–478, 505.
[7]LI W M, SONG P F. A modified ICP algorithm based on dynamic adjustment factor for registration of point cloud and CAD model[J]. Pattern Recognition Letters, 2015, 65: 88–94.
[8]GENG N, MA F F, YANG H J, et al. Neighboring constraintbased pairwise point cloud registration algorithm[J]. Multimedia Tools and Applications, 2016, 75(24): 16763–16780.
[9]程军, 冯丹, 梅魁志. 基于点云配准ICP 算法构建K-D 树的方法: CN110097581B[P]. 2021–02–19.
CHENG Jun, FENG Dan, MEI Kuizhi. Method of constructing K-D tree based on point cloud registration ICP algorithm: CN110097581B[P].2021–02–19.
[10]BESL P J, MCKAY N D. A method for registration of 3-D shapes[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 1992, 14(2): 239–256.
[11]LI Y D, GU P H. Free-form surface inspection techniques state of the art review[J]. Computer-Aided Design, 2004, 36(13): 1395–1417.
[12]薛珊, 吕旸, 吕琼莹, 等. 基于点云几何信息改进的自动配准方法[J]. 制造业自动化, 2020, 42(1): 83–87, 97.
XUE Shan, LÜ Yang, LÜ Qiongying, et al. An improved automatic registration method based on point cloud geometric information[J].Manufacturing Automation, 2020, 42(1): 83–87, 97.