飞机设计性能的不断提升对飞机制造提出了更高的装配要求，出现了基于实测模型的预装配分析技术[1–2]，该技术以实际测量点云数据重构的CAD模型代替理论模型进行虚拟装配分析。模型重构的前提是需要将实测点云划分为只含有单一几何特征的拓扑结构区域。但飞机装配零部件数量多、结构复杂[3]，对大规模点云数据按照设计意图开发智能自动分割算法难度很大，因此如何精确且高效地实现点云分割成为飞机预装配分析研究中亟待解决的问题之一。
国内外对点云分割技术进行了许多研究。董明晓等[4]提出了一种基于数据点曲率变化的区域分割方法，即先对每一条扫描线上的数据点求取曲率值，然后将其中曲率值变化较大的点提取出来作为边界点，将封闭边界包围的区域作为分割的结果，但该方法曲率阈值选取比较困难且只能适用于扫描线格式的点云数据。Besl等[5]最早提出了基于面的区域生长分割方法，该方法从具有特征的某一或多个数据点开始，围绕具有相似几何性质的相邻点生长，重复向外搜索直到分割出所有满足该特征类型和生长参数的点。但需要人工交互选取种子点、确定生长初始参数以及区域生长允许误差。k–means聚类分割算法[6]以空间中k个点为中心进行聚类，对最靠近他们的对象归类。通过迭代的方法，逐次更新各聚类中心的值，直至得到最好的聚类结果。但需要预先估计并给定聚类数量，对初始聚类中心点的选取也较为敏感。赵夫群等[7]提出了基于改进随机抽样一致性 （Random sample consensus，RANSAC）的点云分割算法，通过改进初始点数据的选取方式和判断准则使RANSAC算法对点云数据分割的准确度提高，但仅能从点云数据集中分割出平面数学模型，不适用于含有多种几何特征的复杂点云数据分割。并且针对不同密度信息的点云数据，判断阈值δ需要人工计算后设置。Hao等[8]提出了基于CAD模型驱动先验模板的曲面重构方法，即给定一个CAD模型，提取其中的几何元素，以几何元素为基本单位进行网格划分，构建基于CAD模型的先验模板。并提取网格节点处的法向量，对引导点法线附近的原始扫描点云数据进行搜索，得到的搜索点的高斯加权平均值代表零件在引导点处的实际几何参数。
基于以上分析，本文结合用于虚拟预装配分析的飞机零部件具有理论数模的特点，以具有理论数模的零部件实测点云数据为研究对象，提出了一种基于特征模板的点云数据精确分割技术方法。该方法通过提取理论曲面并构建特征模板实现由特征模板含有的原始信息驱动规模点云数据精确分割出属于目标曲面特征的点云数据，技术路线如图1所示。本文提出的方法已在CATIA的二次开发平台CAA（Component application architecture）上开发了相关程序，并通过后段肋实例验证了方法的有效性和精确性。
图1 基于特征模板的点云精确分割技术路线
Fig.1 Technical route of precise point cloud segmentation based on feature template
1 理论曲面特征模板构建
当前飞机零件的曲面特征多采用自由造型设计，其数学表达式的求解不仅烦琐且效率低，难以将其直接存储在数据库中[9]，因此通过对理论模型待识别曲面特征构建特征模板，将其解构为数据库可存储模式。特征模板构建过程如图2所示，首先对理论曲面特征进行提取命名并放入特征提取集中，对特征提取集中的所有理论曲面特征逐一离散获得离散点集，将离散点作为顶点三角网格化并计算每个三角面片的法向量。最终由曲面离散点集合Pi、三角网格集合Mi、三角面片法向量集合nF共同组成特征模板并存储至特征模板数据库中。
图2 理论曲面特征模板构建过程
Fig.2 Theoretical surface feature template construction process
特征模板的构建过程简单方便且灵活，编号相同的零件能够共用同一组特征模板，减少了特征模板重复构建的工作量，降低了数据库的内存使用率。
1.1 曲面特征离散
曲面特征离散过程如图3所示，先按照设定的离散间距参数采用截交线法[10]将理论曲面特征依次离散为一组截交线，得到的截交线有序集合为L（S） = {Li，i=1，2，…，nL}，对截交线Li通过等弧长法[11]离散，保证相邻离散点之间的弧长相等，并且计算过程简单、效率高。通过CATIA软件开放的二次开发接口函数CATCrvParam获取理论截交线的参数方程，得到的参数表达式为
图3 曲面特征离散
Fig.3 Surface feature discretization
则曲线弧长的计算方程为
给定初始弧长l0，完成单条截交线L的数字化离散，得到有序离散点集合为
整个理论曲面特征离散点集合为所有截交线有序特征离散点的并集，即
式中，i为曲面特征离散截交线的条数；j为第i条截交线上的离散点个数；nP为截交线上总离散点数；nL为曲面上总离散截交线条数。
1.2 三角面片法向量计算
采用三角化算法将曲面特征有序离散点作为三角形顶点连接成三角网格。三角网格是计算机图形学中用于描述不规则实体建立三维数字化模型的一种数据结构，可以表示为顶点、边、面的集合M = [V，E，F]，其中V为网格顶点的集合，E为网格边的集合，F为网格面的集合。
式中，nV为三角网络顶点个数；nE为三角网格边的个数；nF为三角网格面的个数。
为了给后续点云分割提供精确的初始参数以及生长依据，需要提取三角网格模型的面片法向量。三角面片法向量有两个方向，为了避免出现歧义，需要保证三角面片法向量方向的唯一性，如图4所示。
图4 三角面片法向量定义
Fig.4 Triangular face normal vector definition
设两个三角面片分别为f1（P21，P22，P11）和f2（P12，P11，P22），以三角面片f1为例，非公共边的顶点P21作为起始点，同一条截交线的另一个顶点P22作为终点，构建向量
根据右手定则和几何中向量积的定义，由三角形边的向量积求出三角面片f1法向量并单位化。
理论曲面特征的三角面片法向量集合为
2 特征模板驱动的点云数据粗分割
2.1 点云拓扑关系建立
点云数据通常具有大规模和无序的特点，所以建立点云间的拓扑关系对点云数据分割至关重要。根据本研究使用的特征模板，为实现特征模板驱动的点云数据粗分割目标，采用空间栅格法[12]建立点云拓扑关系。空间栅格法是利用空间立体栅格来近似表示点云的空间结构和几何形状，基本原理如图5所示：设xmax、ymax、zmax、xmin、ymin、zmin为点云中x、y、z的最大、最小坐标，给定栅格扩展阈值ε，以点（xmax+ ε，ymax+ ε，zmax+ ε）和点（xmin– ε，ymin– ε，zmin– ε）为对角点且表面平行于空间坐标系的空间六面体构建点云的空间包围盒，将空间包围盒作为栅格划分的总空间，以预设宽度λ作为单位长度并对包围盒进行均等划分，其中栅格划分宽度λ的取值根据数字化扫描设备设置的点分辨率确定，栅格扩展阈值ε可设置为λ的0.5倍。沿x、y、z方向的栅格数量分别为
图5 空间栅格划分原理
Fig.5 Principle of spatial grid division
符号[]表示向上取整。三维立体栅格划分完成后，每个栅格是一个微小正方体；用栅格坐标 （u，v，h）表示所有栅格，u、v、h分别代表在x、y、z方向的栅格序号，则有u∈[0，l – 1]，v∈[0，m – 1]，h∈[0，n – 1]。然后将点云数据点按坐标值放入对应栅格内，有点的栅格称为实格，无点的栅格称为空格；栅格之间的邻域关系间接代表了点云的拓扑关系。
2.2 实格质心点粗分割
为提高点云数据处理效率，并充分利用特征模板的理论信息，将三维点云数据与理论数模对齐后，由特征模板驱动包含有多种特征的原始点云完成粗分割，剔除不属于特征模板的点云数据。具体算法步骤如下：
（1）通过空间栅格法建立全部点云拓扑关系；
（2）按照特征模板中三角面片顶点集合计算并构建曲面特征空间包围盒{Bi（S），i=1，2，…，n}，S为提取的曲面特征，n为特征模板的数量；
（3）建立空的点云数据粗分割块{Ci，i=1，2，…，n}，其中Ci用于存放栅格坐标；
（4）检索空间栅格Gk（u，v，h），若Gk是实格，则以质心点代表栅格实格内所有点，其质心点坐标为
x，y，z）；
（5）遍历所有曲面特征空间包围盒，若栅格的质心点Gk（u，v，h）能成功投影到某一三角面片fi上且投影距离小于εd，则认为质心点
（u，v，h）在曲面特征空间包围盒Bi（S）内部，即
（x，y，z）∈Bi（S），并把栅格坐标Gk（u，v，h）存入对应的分割块Ci中；
（6）重复步骤（4）和（5），直到空间栅格全部遍历完毕，完成点云数据实格质心点粗分割。
曲面特征的空间包围盒具体构建方法如下。
设特征模板三角网格M={V，E，F}中的任一三角面片fi =（Vi1，Vi2，Vi3），Vi1、Vi2、Vi3分别为三角面片fi的3个顶点。将（Vi1，Vi2，Vi3）分别沿nF（fi）和–nF（fi）方向移动εd构成正偏距三角面片fi′=（Vi1′，Vi2′，Vi3′）和负偏距三角面片fi″=（Vi1″，Vi2″，Vi3″），其中nF（fi）为三角面片fi的法向量，εd为包围盒容忍阈值。以fi′和fi″为上底面和下底面，Vi1′Vi1″、Vi2′Vi2″、Vi3′Vi3″为3条棱构成的封闭三棱柱作为单个三角面片fi的空间包围盒B（fi），即
曲面特征的空间包围盒B（S）为特征模板中所有三角面片包围盒的并集组成，即
式中，n为三角面片的个数。
如图6所示，
分别为不同空间栅格实格的质心点，其中认为
在空间包围盒Bi（S）内部，而因为
不能成功投影到三角面片fi上、
在三角面片fi投影距离超过阈值εd认为
均在空间包围盒Bi（S）外部。
图6 空间包围盒建立与质心点粗分割
Fig.6 Spatial bounding box establishment and coarse division of centroid points
3 基于面的自适应区域生长细分割
对获得的粗分割点云块进一步进行细分割，由特征模板数据库离散点集合确定理想种子点，理想种子点的理论法向量作为生长分割初始参数向外扩展，搜索种子点的k近邻点，采用主成分分析方法对点云中每一点法向量进行估计，并基于邻域方向一致性传播的调整方法来统一所有点云法向估计的方向，若满足过渡夹角阈值条件的点作为生长点继续生长，若生长到点云边界则停止生长。
3.1 种子点选取
读取存储在数据库的特征模板对应的离散点集信息，选择Pi=P[nL/2][np/2]作为细分割理想种子点Pseed。其中nL为曲面上总离散截交线条数；np为截交线上总离散点数；[]表示向下取整。Pseed为有序离散点集的中间点，保障了其所属点云区域的准确性，确保种子点从曲面特征中心由内向外开始生长。
同时根据特征模板数据库索引关系检索理想种子点Pseed作为共享顶点的所有三角面片的编号，所有三角面片法向量的矢量和单位化后作为Pseed的法向量nseed（图7），赋予Pseed精确的初始生长分割参数，能够有效利用理论曲面特征信息引导点云精确分割。
图7 理想种子点法向量定义
Fig.7 Definition of ideal seed point normal vector
3.2 点云法向量估计
为了给后续点云数据细分割提供生长依据，需要对点云法向量进行估计。本研究采用基于局部表面拟合的方法[13]对大规模散乱点云的法向量进行初始估计。
对于点云{Vi，i=1，2，…，n}中的某一点xi∈{Vi}，通过空间栅格法建立的拓扑关系搜索点云中距离点xi最近的k个点，即点xi的k邻域，记作kNb（xi），对kNb（xi）最小二乘拟合局部切平面P，P的法向量记作n，点xi到坐标系原点的欧式距离记作di，则局部切平面P可以表示为
采用主成分分析方法 （Principal component analysis，PCA）[14]求解切平面法向。对式（13）的求解问题就可以转化为对式（14）中的半正定协方差矩阵M的特征值求解问题。
式中，
为kNb（xi）的质心点，将M分解为奇异值，得到的最小主成分系数作为P的法向量n，使得
式中，U为左奇异矩阵，由正交向量组成；Σ对角线的元素为奇异值；V为右奇异矩阵，由正交向量组成。V的列向量是变换后的新的主成分向量，即两个相切向量和一个垂直于曲面的向量。
通过主成分分析法计算得到的点云法向量有可能指向点云模型外部或者内部，为了避免产生歧义，需要对点云法向量的方向进行一致性调整。本文采用基于邻域方向一致性传播的调整方法：首先以点云中种子点Pseed的法向量nseed为基准，点Pseed的k邻域为kNb（pseed）={pi，j=1，…，k}，点Pseed的相邻点pi的单位法向量为ni，两点之间的法向量的夹角为θ = arccos（nseed·ni），0≤θ≤π，θ越小表示两点的法向量方向一致程度越高。设置法向量夹角阈值θmax，若θ∈[0，θmax]，表示两点法向量方向是一致的，则ni保持不变；若θ∈[π – θmax，则把ni反向，得到调整后的法向量ni；若θ∈（θmax，π – θmax），表示两点法向量方向没有关联。种子点Pseed的k邻域法向量方向全部调整完毕后，再以pi调整后的法向量为基准对pi的k邻域kNb（pi）中未调整法向量方向的点进行调整，不停传播下去直到所有点云数据的法向量方向调整完毕。具体原理如图8所示。
图8 点云法向量一致性调整
Fig.8 Point cloud normal vector consistency adjustment
图9为采用主成分分析法对Fandisk点云模型进行法向量估计的结果，从局部放大图可以看出，点云在光滑区域、边界处和拐角处的法向量估计结果均接近理论状态，且经过法向一致性调整能实现点云法向量均指向Fandisk点云模型外侧。
图9 Fandisk点云模型法向量估计结果
Fig.9 Fandisk point cloud model normal vector estimation results
3.3 自适应区域生长
为了得到精确的点云生长分割结果，避免区域生长分割过程中人工确定初始种子点和反复尝试分割参数，同时结合曲面特征内部局部点云法向量变化相对平滑和边界处点云法向量突变明显的特点，本文提出了基于面的自适应区域生长方法实现点云的细分割，即通过栅格实格质心点粗分割过程实现与其所属三角面片及法向量建立索引关系，便于后续实格质心点归属关系查询。从选定的理想种子点Pseed开始，若相邻两点的法向量夹角小于各自归属的三角网格面片间的法向量夹角，则作为新的种子点继续生长；伴随着生长点不断向外搜索扩张，允许生长的过渡夹角阈值参数也随着曲面三角面片的几何性质和类型不断自适应调整，无需通过已生长点逼近不同类型曲面识别特征属性；若生长点的法向量与所属的三角面片法向量夹角大于投影夹角阈值则认为区域生长到达了边界处，如图10所示，曲面S点云内部的相邻两点P1和P2过渡夹角较小，但边界点P3的投影夹角较大。
图10 曲面特征点云过渡夹角和投影夹角分布
Fig.10 Distribution of transition angle and projection angle of feature point cloud on curved surface
基于面的自适应区域生长点云细分割曲面特征细分割具体步骤如下。
定义种子点容器Seed和细分割点云块集合C细。
（1）原始点云经过粗分割得到的粗点云块C粗和选取细分割理想种子点Pseed，在C粗中找到距离Pseed最近点作为初始种子点p1并加入种子点容器Seed中。
（2）从Seed中的种子点p1开始，通过空间栅格法建立的点云拓扑关系查找点p1的k近邻点pk，并遍历所有k近邻点pk。
（3）计算pk内的点与点p1的法向量偏差en，其中过渡夹角en = arccos（np1·npk），np1和npk分别为点p1和pk估计的单位法向量。
（4）计算pk内的点法向量与其归属的三角面片的法向量nF（fi） 偏差ep，其中投影夹角ep = arccos（npk·nppk），npk和nppk分别为点pk的单位法向量和其归属三角面片法向量。
（5）若pk C细且同时满足en < ε（en）和ep < ε（ep），则把pk同时放入C细和Seed的末位；若pk C细但只满足en < ε（en），则只将pk放入C细的末位。ε（en）和ε（ep）分别为过渡夹角阈值和投影夹角阈值。
（6）遍历完种子点p1的所有k邻域点pk后，从Seed中删除p1。
（7）重复执行步骤（2）～（6），直到Seed内无种子点，点云细分割完毕。
4 软件系统开发与实例验证
在航空工业中，国内大多数企业采用CATIA三维数字化软件，因此本文在CATIA二次开发平台CAA中开发了基于特征模板的点云数据精确分割软件并实现了预期功能，软件界面如图11所示，共包含曲面特征提取模块、特征模板构建存储模块和点云精确分割模块3个模块。
图11 基于特征模板的点云精确分割软件界面
Fig.11 Software interface for precise point cloud segmentation based on feature template
试验测试中，使用激光跟踪仪+T–Probe测量后段肋零件工装定位面上不共线的3个点建立测量坐标系OM，激光跟踪仪+T–Scan扫描后段肋零件表面获取的点云数据位于测量坐标系OM下，将实测点云数据xi转换到理论CAD模型所在零件坐标系OP下以实现与理论CAD模型的对齐。
式中，xi'为转换后的点云数据；R和T分别为测量坐标系OM变换到零件坐标系OP的旋转矩阵和平移向量。
图12给出了应用本研究点云数据精确分割软件对某飞机后段肋实测点云的精确分割实例。后段肋零件理论三维模型如图12（a）所示；图12（b）对后段肋所有待识别曲面特征进行提取命名并加入到特征提取集内；图12（c）将特征提取集内的所有曲面特征构建成特征模板并将构建过程获得的曲面离散点集合、三角网格集合和三角面片法向量集合等信息储存至特征模板数据库中；图12（d）显示了测量坐标系和零件坐标系在空间中的相对位置关系；图12（e）为坐标系对齐后实测点云数据的偏差分析彩图，所含点数85456；图12（f）显示了划分空间栅格的参数设置和建立点云间拓扑关系的具体过程；由特征模板驱动的点云粗分割结果如图12（g）所示；图12（h）给出了粗分割点云法向量估算和一致性调整后的结果；图12（i）给出了实测点云数据精确分割的最终结果，用不同颜色区分识别不同曲面特征的点云数据并注明每个点云分割块内所含点数。结果表明特征模板内的所有待识别特征点云全部完成精确分割。
图12 后段肋实测点云数据精确分割实例
Fig.12 An example of accurate segmentation of measured point cloud data of posterior rib
5 结论
本文提出了一种基于特征模板的点云数据精确分割方法。用该方法对某型飞机后段肋零件实测点云中对待识别曲面特征的点云进行了精确分割。实例验证表明该方法可以按照操作人员的意愿定制特征模板，并通过特征模板驱动和引导分割出点云数据中属于目标曲面特征的点云数据，在满足分割效率和自动化要求的情况下提高了点云分割准确性，能够得到精确的分割结果。
使用本方法时需要注意：分割过程中依赖点云数据与理论曲面特征模板的空间位置关系判定，因此实测点云数据分割前应与理论三维模型进行配准，并且配准的精度会在一定程度上影响后续点云分割结果。
[1] 包宗明.基于实测模型的舱段虚拟装配与容差优化研究[D].沈阳: 东北大学, 2012.
BAO Zongming.The research on virtual assembly and capacity difference of cabin based on measured model[D].Shenyang: Northeastern University, 2012.
[2] 周石恩.基于数字孪生的复杂产品装配建模与精度分析方法[D].杭州: 浙江大学, 2019.
ZHOU Shien.Assembly modeling and accuracy analysis method of complex product based on digital twin[D].Hangzhou: Zhejiang University, 2019.
[3] 薛翔, 张彤阳.现代飞机装配工艺及展望[J].现代制造技术与装备, 2018(10): 184–185.
XUE Xiang, ZHANG Tongyong.Modern aircraft assembly process and prospect[J].Modern Manufacturing Technology and Equipment,2018(10): 184–185.
[4] 董明晓, 郑康平, 姚斌.曲面重构中点云数据的区域分割研究[J].中国图象图形学报, 2005, 10(5): 575–578.
DONG Mingxiao, ZHENG Kangping, YAO Bin.Research on point cloud data segmentation in surface reconstruction[J].Journal of Image and Graphics, 2005, 10(5): 575–578.
[5] BESL P J, JAIN R C.Segmentation through variable-order surface fitting[J].IEEE Transactions on Pattern Analysis and Machine Intelligence, 1988, 10(2): 167–192.
[6] HARTIGAN J A, WONG M A.Algorithm AS 136: A K–means clustering algorithm[J].Applied Statistics, 1979, 28(1): 100–108.
[7] 赵夫群, 马玉, 戴翀.基于改进随机抽样一致的点云分割算法[J].科学技术与工程, 2021, 21(22): 9455–9460.
ZHAO Fuqun, MA Yu, DAI Chong.Point cloud segmentation algorithm based on improved random sample consensus[J].Science Technology and Engineering, 2021, 21(22): 9455–9460.
[8] HAO L, HUANG X, LI S G.Surface reconstruction based on CAD model driven priori templates[J].The Review of Scientific Instruments, 2019, 90(12): 125116.
[9] 毛喆.飞机装配几何特征三维测量模型构建技术研究[D].南京: 南京航空航天大学, 2020.
MAO Zhe.Building technology of 3D measurement model for geometric features in aircraft assembly[D].Nanjing: Nanjing University of Aeronautics and Astronautics, 2020.
[10] YA 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.
[11] 刘丹丹, 张树有, 刘元开, 等.一种基于特征点识别的曲线离散化方法[J].中国图象图形学报, 2004, 9(6): 755–759.
LIU Dandan, ZHANG Shuyou, LIU Yuankai, et al.A curve discretization method based on recognization of feature point[J].Journal of Image and Graphics, 2004, 9(6): 755–759.
[12] 张春成.一种改进的点云数据栅格算法[J].兵工自动化,2013, 32(6): 90–93.
ZHANG Chuncheng.An improved grid algorithm of point cloud data[J].Ordnance Industry Automation, 2013, 32(6): 90–93.
[13] AMENTA N, BERN M.Surface reconstruction by voronoi filtering[J].Discrete & Computational Geometry, 1999, 22(4): 481–504.
[14] ABDI H, WILLIAMS L J.Principal component analysis[J].Wiley Interdisciplinary Reviews: Computational Statistics, 2010, 2(4):433–459.