1 引 言

激光扫描技术分为地面TLS和机载LiDAR两种。机载LiDAR在数字城市方面得到了广泛应用，可用于提取城区的树木^[1]、城区地物的分类及特征提取^[2]。TLS主要用于航空、航天、船舶、汽车和模具等现代制造业邻域，也可应用于修复破损的艺术品或者被损零件^[3]，例如修复破损的雕像、雕刻；修复由于常年使用造成的船舶破损；修缮古建筑等。此时并不需要对整个原型进行复制，而只需借助于逆向工程技术提取需要修复的零部件以实现缺损空洞的定位和逆向建模，仿真缺损物体的修理工艺，实现修复技术的现代化。在研究修复过程中，需要对三维激光扫描测量技术获取的点云数据进行空洞修复研究。在实际应用中空洞有不同的种类，按曲率大小可分为大曲率空洞和小曲率空洞；按空洞大小可分为大空洞和小空洞；此外还可分为封闭空洞和非封闭空洞，以及有特征和无特征空洞。

目前空洞边界自动提取方法主要有两种：第1种是在曲面网格模型已经构建的基础上，通过网格的拓扑关系识别出网格中存在的空洞。这种方法的主要思想是：如果网格模型中没有空洞，那么模型中的任意一条边都属于且仅属于2个三角片，如果存在某条边仅属于1个三角片，则该条边是构成空洞多边形的1条边，即为边界边^{[4, 5]}。根据边界边在网格模型中的定义，只需遍历整个三角网格面就可以找出空洞边界。文献[6]按照三角片的边长情况，对用边界连接得到的空洞三角片细分，并采用类似Delaunay三角化的思想，调整相邻三角片的公共边，在空洞内部不断生成新增三角片。文献[7]则采用体数据场融合的方法进行空洞修补。该方法比较简单，但其前提是需要预先生成三角网格模型，而在三角网格模型的生成过程中有可能已经忽略了空洞特征的存在^[8]。第2种是直接从散乱点云中识别出边界点。文献[9]提出首先计算每个点的邻近点，确定空洞边界，然后利用边界点的邻近点构造曲面来填补空洞。文献[10]根据空洞周围的局部离散点来建立一张曲面片，然后采用面上取点的策略补出空洞部位所缺的点。文献[11]提出一种通过判断k近邻域分布均匀性识别空洞边界，其实是利用了数据挖掘方法中的聚类分析通过角度来判断边界点^[12]，该方法由于噪声的干扰会使判断出现错误；邱泽阳提出了一种将散乱点云投影至二维平面上提取边界点的方法^[13]，由于会出现投影重叠等问题，所以不能很好地应用于大曲面修复中。

基于以上原因，本文设计了一种推进式逐层求解法ALS提取空洞边界，首先将求解出的k近邻域利用ALS投影至平面，然后在平面内进行网格划分，提取边界网格，在此基础上应用最小凸包法提取边界线；最后，进行空洞边界与物体本身边界的识别。

空洞边界提取需要利用空洞附近数据点的拓扑关系以及这些数据点在平面的投影点，然后提取边界网格，具体步骤如下：

(1) 搜索k近邻域。这是逆向工程中点云数据预处理中常见问题之一，本文采用改进式的非均匀栅格法构建空间点云数据，缩小搜索范围，从而节约计算时间，提高算法效率。

(2) 推进式逐层求解投影。为了保留特征和避免投影重叠，本文将计算出的近邻域点采用推进式逐层求解投影方法。

(3) 提取边界网格。利用投影点进行网格划分，根据实网格和空网格邻近关系的判定，提取空洞边界网格。

(4) 最小凸包法提取边界线。针对每一个网格，应用最小凸包法提取边界线，避免噪声点的干扰。

(5) 边界线识别：根据提出的边界线进行识别，即将需修复的空洞边界与物体本身存在的边界进行区分。

具体提取流程如图 1所示。

图 1 空洞边界提取流程图 Fig. 1 Flowchart of hole boundary extraction

图选项

本方法以没有任何拓扑信息的散乱点云为研究对象，因而首先需要建立点的拓扑关系。点云的拓扑关系定义为k近邻域，某点的k个近邻域的计算是指在数据集S={p_i，i = 1，2，…，n}中找到k个与该点欧氏距离最近的点^{[14, 15, 16, 17]}。为此，许多学者对此进行了一些快速算法的研究，如k-d tree法^[18]、八叉树算法^[19]。这些方法对于在大曲面修复中分布不均匀的海量数据而言显得效率较低。为此本文采用改进式的非均匀栅格法，空间数据点集总数N，k近邻域，数据点的密度ρ等因素综合制约着栅格分割及k近邻域的搜索。具体算法步骤如下：

(1) 根据最小包围盒的思想，利用式(1)构造被测曲面的最小外接立方体；

点击浏览公式

点击浏览公式

点击浏览公式

(3) 对所划分的小立方体进行3个坐标方向编号，从而确定在小立方体内数据点所属的编号；

(4) 根据已知点P(x,y,z)，求得其对应子立方体编号以及相邻的26个子立方体编号；

点击浏览公式

(5) 在计算出的27个子立方体内部搜索，找到距离最近的k个近邻点，如果点数不够，则将子立方体向外扩展，继续搜索，若仍达不到要求，则将P点认定为噪声点删除，继续搜索下一个点，直到所有的点搜索完毕，得到每一个点的k近邻域。流程具体如图 2所示。

图 2 k近邻域搜索算法流程图 Fig. 2 Flowchart of k neighborhood search algorithm

图选项

推进式逐层求解法是利用空间数据的多尺度特征，此特征是指空间数据在不同观察层次上所遵循的规律以及体现出的特征不尽相同^[20]。此方法的主要思想是利用同一层次上的点云数据具有相同曲率这一特征而提出的推进式逐层求解法。通过逐层求解各点的k近邻域点的曲率。曲率可以看作是法矢沿给定方向的变化率^[21]，法矢的计算可参考文献[22] ，假设其中某一点p的法矢为N，指定一方向矢量(Q-P)为从点p指向它邻域上的一点q，则在点p处沿该方向的曲率为：

点击浏览公式

最小曲率k_min(p)=k( p,q₁)=

最小曲率方向D₁ = Q₁-p

最大曲率

最大曲率方向D₂ = Q₂-p

平均曲率

根据平均曲率设定1个阈值γ，把曲率大于阈值的点划分为1层，然后对同一层的点计算其中点c，为了体现c的影响。为此给任意点p_m赋予1个权值，这个权函数为高斯函数

点击浏览公式

点击浏览公式

点击浏览公式

点击浏览公式

求得投影平面后即可进行三维数据到二维的投影。根据法矢夹角准则确定邻域投影点，判断点法矢n_p与邻域点法矢n_i的夹角θ，当θ < θ_t，θ_t为夹角阈值，一般设θ_t≤π/2，此时当三维的点投影到二维平面上时，位置关系不发生改变，可以提高投影点数；根据得到的交点数据集进行坐标投影，将三维坐标点投影到平面上。给定散乱数据点集，赋予每个点投影重复度属性，可以用POINT.project表示，初始值设为零，被投影一次则POINT.project++，这样在选择点p投影时，首先对其重复度属性判断，如有重叠投影，则进行重新分层。将所有点的k近邻域采用以上方法进行计算，其中某两点得出的投影效果如图 3所示。图 3(a)为没有空洞边界时k近邻域的投影图，图 3(b)为有空洞时的投影图。

图 3 k近邻域点投影图 Fig. 3 k nearest neighbor projection

图选项

在投影平面内，找到整个点的k近邻域的最小包围盒，用一定间隔(由k值决定)矩形网格将包围盒剖分，如图 4所示。该网格的网孔可按其是否拥有散乱点分为2类：一类为实孔；另一类为空孔。每个网孔有4个相邻的网孔(除网格边界外)，对于实孔，若它的4个相邻网孔中至少1个为空孔或它的相邻网孔少于4个时，认为它是边界孔(图 4中的黑网格)，将提取出的所有边界孔进行连接，就可得到边界网格。

图 4 提取出的边界网格 Fig. 4 Extracted boundary grid

图选项

进一步的工作就是利用最小凸包法在边界网格内连接出空洞边界曲线，具体处理方法如下：

(1) 4个相邻网孔中只有1个空孔的情形，将与空孔相邻的边称为空边。若空边与y轴平行，则将边界网格内散乱点中的y_max与y_min的点连接成直线AB，将与其空边的一侧散乱点作最小凸包，除去AB的凸包边界即为空洞边界，空边平行于X轴的处理方法类似，详见图 5。

图 5 一条空边时的空洞边界提取 Fig. 5 The hole boundary extraction with an empty side

图选项

(2) 4个相邻网孔中有2个空孔的情形，第一条边与只有一条空边的情形有相同处理，在另一区域内处理第二条边；对于具有2条不相连的空边情形，先将区域分为2个部分，然后依据上述方法进行处理，处理结果如图 6所示。

图 6 2条空边时的空洞边界提取 Fig. 6 The hole boundary extraction with two empty sides

图选项

(3) 4个相邻网孔中有3个空孔的情形，求1个网格内所有散乱点的最小凸包，将靠近实边的2点断开，其余边则为空洞边界，如图 7所示。

图 7 3条空边时的空洞边界提取 Fig. 7 The hole boundary extraction with three empty sides

图选项

将所有的网格按上述方法进行处理，连接所提取的边界得到整个边界，而这个边界是空洞边界还是物体本身的边界，还必须进行识别。识别的原则是：以曲面外法矢方向为Z轴正方向，若规定沿边界边遍历时其邻域点在它的右侧，则如果边界线的遍历方向是逆时针，则此边界线为空洞边界；反之此边界线为物体本身的外边界。

为了验证本文方法的稳健性和可靠性，将零件模型挖出1个空洞进行实验。实验数据为真实物体的三维点云数据，图 8表示用文献[11]和本文推进式逐层求解方法提取的边界线。从图中比较可知，文献[11]对于这种二值曲面只能识别出存在1条空洞边界线。为了更能体现本文方法的优越性，将所提取的空洞边界线应用于后续的空洞填充和曲面重构中，图 9为文献[11]提取的空洞边界应用效果图，从图中可以看出模型的特征信息不能很好地被保留，而且反面模型仍然存在空洞。

图 8 带空洞的零件模型 Fig. 8 The part model with a hole

图选项

图 9 文献[11]提取的空洞边界的应用 Fig. 9 The application of extracted hole boundary for literature[11]

图选项

图 10为本文推进式逐层求解法提取的空洞边界应用效果图，由于该方法是基于曲率进行逐层求解的，所以物体的特征信息能很好地保留下来，且正反面模型都修复完好。通过比较证明本文提出的方法具有很好的应用效果，能满足后续工作的需要。为了进一步验证本文的方法，将该方法应用到某曲面模型的修复，图 11是将本文的方法进行边界的提取和填充的效果图。图 12表示了物体内外边界的识别，图中箭头方向表示在识别边界时的走向，粗线表示物体本身的边界而细线表示空洞边界。

图 10 本文方法提取的空洞边界的应用 Fig. 10 The application of extracted hole boundary for this method

图选项

图 11 本文方法提取的空洞边界及填充 Fig. 11 The application of hole boundary and filling for this method

图选项

图 12 内外边界识别 Fig. 12 Identification of inside and outside boundary

图选项

试验中，对物体分别设置的k、α值的组合进行试验，结果如表 1所示，可以看出不同数量的点云数据和不同近邻数k的最快搜索速度的α值范围比较集中，一般在1～2之间，当k值较小时，α值应取大些；当k值较大时，α值应取小些，这样可控制子立方体边长过大或过小的情形。

为了能定量地评价空洞边界提取的精确性，本文将挖出的空洞点云作为采样控制点，计算控制点至修补曲面模型的距离，用最大误差和平均误差作为精度评价指标。本文采用空间快速迭代算法^[23]来计算点到曲面的最小距离。如图 13所示，Q为修补曲面S(u,v)外一点，将曲面进行网格划分，与Q距离最近的网格点为初始点P₀，P_c是Q在曲面上的投影点，是点Q到曲面的最小距离，记为d，n是曲面在P_c处的法矢量。

图 13 点到曲面的距离 Fig. 13 The distance from point to surface

图选项

于是可得

点击浏览公式

用S_u、S_v分别乘式(11)两边，由于n=,所以S_und=S_vnd=0，设E=S_uS_u F=S_u·S_v，G=S_v·S_v,可得

点击浏览公式

解上述方程组，可得Δu、Δv，令u+Δu→u，v+Δv→v，重复以上过程，直到|Δu| < ε，|Δv|<ε，此时，即得点Q到曲面的距离d。空间快速迭代算法效率很高，一般迭代3~10次即可收敛，而且在边界处的稳定性很好，可以满足需要。

得到点到曲面的距离后，遍历所有采样点，可得到1个集合{d_i}₁ⁿ，n为采样点数，将max{d_i}ⁿ₁定义为采样数据到曲面模型的最大误差;将定义为采样数据到修补曲面模型的平均误差。计算结果如表 2所示。

通过对本文推进式逐层边界提取方法与传统边界提取方法对比，证明本文提出的方法在噪声干扰和多值曲面条件下具有良好的稳健性和适应性。研究结果表明：① 本文通过曲率的计算，在提取边界时保留了其特征信息，所以在该边界约束下填充的点云能更好地拟合缺失数据；② 本文方法利用数据的多尺度特征进行分层，避免了投影重叠的问题，从而可以适用于多值曲面的修复；③ 本文在边界网格内用最小凸包法精确提取边界线，能很好地避免噪声干扰，使提取的边界更精确。因此，针对多值曲面的简单空洞具有适用性。由于所采集的数据不同，面临的空洞曲面类型也不同，处理的方法也不同，本文对空洞区域存在多种曲面类型的情况仍需进一步探究，这是后续研究的重点。