基于线激光的航空叶片孔洞修复方法

A Method for Repairing Holes in Aviation Blades Based on Line Laser

摘要: 针对航空叶片的孔洞修复问题，本文提出了基于线激光的孔洞修复方法。使用Alpha Shape边缘轮廓提取算法提取航空叶片的内外边缘，使用基于KD树的欧几里德分割算法分割出孔洞部位的边缘，接着使用基于最小二乘的点云曲线拟合算法，根据线激光扫描的点云线条拟合曲线方程，再将孔洞区域的边缘数据导入拟合后的曲线模型得到孔洞部位曲面数据，完成航空叶片孔洞的修复工作。

Abstract: For the repair of holes in aviation blades, this article proposes a hole repair method based on line laser. Using the Alpha Shape edge contour extraction algorithm to extract the inner and outer edges of aviation blades, using the KD tree based Euclidean segmentation algorithm to segment the edges of the hole area. Then, using the least squares based point cloud curve fitting algorithm, fitting the curve equation based on the point cloud lines scanned by line laser, and importing the edge data of the hole area into the fitted curve model to obtain the surface data of the hole area, Complete the repair work of aviation blade holes.

1. 引言

航空发动机叶片是发动机中最关键的高性能部件之一 [1] ，其外形精度和表面完整性控制着发动机的性能水平和使用寿命。航空发动机叶片在与发动机主轴、涡轮盘和腔体配合的过程中，改变空气或燃气的方向，形成高速、高温的易燃压缩气流，保证飞机所需的动力。由于叶片长期在高温、高压的恶劣工作环境下高速运转，叶片的损伤在所难免，约70%的发动机故障是由叶片故障引起的。因此，航空叶片的修复再制造显得至关重要，而对叶片的孔洞修复则是这个环节中经常需要面临的难题。

对于孔洞的修复，谢文军等 [2] 针对构建实景的孔洞问题，使用边界收缩以及构建三角面片的方法对孔洞进行填补。肖俊等 [3] 针对岩体的激光点云的孔洞问题，提出了构建孔洞坐标系，对孔洞进行新增点采样，再进行二维投影插值修复的方法，熊威等 [4] 针对植物叶片点云空间重叠大且离群点多的问题，使用改进的区域生长算法去除离群点对点云进行三角剖分完成孔洞修复。

然而上述方法在对航空叶片的使用过程中，存在计算量大且无法保证叶片修复后的孔洞与叶片本身弯曲度的一致性。本文根据上述存在的问题，提出了基于线激光的航空叶片孔洞修复方法，该方法能够在较快时间内修复点云，且修复效果满足后续加工的需求。

2. 线激光点云扫描及点云滤波

本文通过线激光轮廓仪进行点云扫描，采用交叠式扫描方法获取涡轮叶片曲面点云，通过点云拼接完成点云点云的扫描流程。扫描后的点云由于有孔洞的存在，采用统计滤波的方式对点云进行降噪。

2.1. 交叠式扫描方法

如图1为线激光扫描系统的扫描流程，红框和绿框为系统前后两次扫描的区域，黄色箭头为线激光轮廓仪的运动轨迹，系统可以控制扫描次数，从右到左按照黄色箭头的运动轨迹，重复n次上下扫描完成叶片的扫描，这样就可以得到叶片的这一面表面轮廓数据。点云数据用txt文件格式保存，方便后面数据的读取，由于数据每次都是以当次上下扫描的起始位置作为坐标原点进行保存，所以n次上下扫描后会后n组不同坐标系的数据，为了将所有坐标系下的数据统一到同一坐标系下，需要对数据按照激光轮廓仪的位移量进行平移。假设第二次上下扫描时向做移动了20 mm，则将第二次扫描的数据向右平移20 mm，以此完成坐标系的统一。但由于系统运动中会有误差，所以为了提高数据精度，每次对平移的相互两幅点云进行ICP配准来减小因误差而产生的间隙。为了ICP配准的精度，需要在前后两幅点云间设置交叠区域，以交叠区域作为特征进行配准，经过多次实验发现，当重叠率大于50%时，点云的配准成功率和配准精度会提高，但是会增加数据采集时间以及配准时间，而当重叠率小于甚至是小于20%时，会产生较大的配准误差以及匹配错误的情况，所以交叠区域的重叠率未50%最为合适 [5] 。

Figure 1. Line laser system scanning process

图1. 线激光系统扫描流程

2.2. 统计滤波降噪

在使用线激光扫描获得叶片点云后，由于孔洞的存在而导致噪声的产生，需要使用统计滤波的方法将噪声滤除。

统计滤波 [6] 是根据点云的邻域信息建立异常点剔除准则。异常值通常是无组织的，远离主点云，稀疏，几何上不连续，局部点密度不一致。相比之下，主点云则相对集中和密集。统计滤波算法是基于离群点与相邻点之间的距离较大，而主点与相邻点之间的距离较小的特点。对每个点邻域的统计分析用于去除异常值。首先，建立一个八叉树，以点云的最小封闭立方体作为根节点，并将该立方体分成8个大小相等的子立方体。如果一个子立方体包含点，那么它应该继续被划分为8个大小相等的立方体，直到点云中的每个点获得一个唯一的索引坐标。然后，计算每个小立方体中点云的密度，另一方面，而且通常，包含异常值的子立方体的密度相对较小。具体的实施步骤如下：

1) 建立一个自适应八叉树栅格化点云和计算的密度云 ${\rho }_n$ 在每个网格。

2) 对于在每个网格中的点云P的任意点 $p_j\in P$ ，搜索 $p_i$ 的k近邻点 $p_i\left(j\in k\right)$ ，并计算 $p_i$ 和 $p_j$ 之间的距离。

3) 将一个点是离群值的概率定义为：

$\delta =\frac{{\displaystyle {\sum }_{j=1}^k\Vert p_i-p_j\Vert }}{k{p}_n}$ (1)

4) 建立一个阈值。如果 $\delta$ 大于阈值，则认为 $p_i$ 是一个离群值并被删除。

3. 涡轮叶片孔洞区域边缘提取

将曲面点云投影到二维平面变成平面点云，使用Alpha Shape算法 [7] 实现边界提取，再将获得的数据还原到三维曲面上，得到曲面点云的边缘轮廓。

3.1. 基于Alpha Shape的曲面边缘提取

本文提出的曲面轮廓提取算法是通过将曲面点云投影到平面在还原的方式实现。如图2(a)为Alpha Shape提取轮廓点示意图，在平面点云中，有一个圆在点云的任意两点之间滚动，这个圆的半径取值为α，圆在点云上滚动时需要满足一定的条件；而当α的取值很小时，点云中的每个点都会被判断为边界点，而当α的值达到一个足够大的阈值时，这个圆只能在平面点云的边界上移动，此时圆移动的轨迹可以被视为边界点。

(a) (b)

Figure 2. Schematic diagram of Alpha Shape algorithm. (a) Diagram of extracting contour points from Alpha Shape; (b) determination diagram of Alpha Shape contour points

图2. Alpha Shape算法提取轮廓点示意图。(a) Alpha Shape提取轮廓点示意图；(b) Alpha Shape轮廓点判断图

其具体的算法步骤如下：

(1) 对于任一点 $p\left(x,y\right)$ ，滚动半径为α的圆，在点云内部搜索到距离p点固定阈值内的所有点，阈值取值为2α，记为点集Q。

(2) 选取点集Q中任意一点，根据和半径的值，计算出两个圆心坐标 $o_1\left(x_{o1},y_{o1}\right)$ 和 $o_2\left(x_{o2},y_{o2}\right)$ ，如图2(b)所示。其坐标计算公式如下所示：

$x_{01}=x+\frac{1}{2}\left(x_1-x\right)-H*\left(y_1-y\right)$ (2)

$y_{01}=y+\frac{1}{2}\left(y_1-y\right)-H*\left(x-x_1\right)$ (3)

$x_{02}=x+\frac{1}{2}\left(x_2-x\right)-H*\left(y_2-y\right)$ (4)

$y_{o2}=y+\frac{1}{2}\left(x_2-x\right)-H*\left(y_2-y\right)$ (5)

其中：

$H=\sqrt{\frac{a^2}{S^2}-\frac{1}{4}}$ (6)

$S^2={\left(x-x_1\right)}^2{\left(y-y_1\right)}^2$ (7)

(3) 在Q点集中，计算除 $p_1$ 点以外的点到 $o_1$ 的距离 $L_1$ 和到 $o_2$ 的距离 $L_2$ 。终止条件为 $L_1$ 或 $L_2$ 的值大于α，则p点为点云的边缘点。

(4) 若 $L_1$ 和 $L_2$ 不全都大于α，则选取另外一点作为 $p_1$ 点。在计算新的圆心坐标之后，重复(3)中的步骤的判断。如果遍历完所有Q邻域的近邻点，有不能满足的(3)的条件，则p点不是点云的边缘点。

本文将曲面点云投影到平面坐标系上，将曲面点云变为平面点云，再使用Alpha Shape算法得到平面点云的边缘。然后获取平面边缘的y，z的坐标，将边缘点云的y，z坐标的值与原曲面点云的y，z坐标进行比对，如果y，z值相同，则认为该曲面上的点为平面边缘点的对应点，从而得到叶片的曲面边缘。

3.2. 欧几里得聚类

提取出内外边缘后，通过点云可视化可以发现，点云的内外边缘的欧式距离很大，所以选择使用欧几里德分割算法来进行点云的内外边缘分割。

欧氏聚类是一种基于欧氏距离度量的聚类算法 [8] 。三维空间中两点之间的欧几里德距离可以定义为：

$d_E=\sqrt{{\left(x_1-x_2\right)}^2+{\left(y_1-y_2\right)}^2+{\left(z_1-z_2\right)}^2}$ (8)

基于KD树的最近邻查询算法是加速欧式聚类过程的前提。对于给定的查询数据点P，需要从KD树的根节点进行比较，其中P(k)是当前节点分区维度k中的数据点P对应的值，m是当前节点分区的阈值。如果为P(k) < m，则访问左侧子树，否则访问右侧子树，直到到达叶节点Q。Q是当前的最近邻，P和Q之间的距离是当前的最小距离。然后将原始搜索路径追溯到根节点。如果在此过程中发现一个距离小于D的点，则需要包含搜索类别中未访问过的所有子节点，并及时更新最近的点，直到所有搜索路径都为空为止。此时，基于KD树结构的整个最近邻查询过程将被声明为已完成。对于欧氏聚类，距离判断准则是前面提到的欧氏距离。对于空间中的一个点p，通过KD树k个最近邻搜索算法找到最接近点p的k个点，以及距离小于集合阈值的点被聚类到集合q中。如果Q中的元素数量不再增加，则整个集群过程将结束。否则，就需要在集合Q中选择除P点以外的点，并重复上述过程，直到Q中的元素数量不再增加为止。

4. 基于最小二乘的涡轮叶片点云曲线拟合

因为多项式曲线有便于求导的优势，又可以写为参数方程的形式易于采样，本文使用最小二乘法 [9] 多项式拟合孔洞线。使用最小二乘进行曲线拟合就是寻找一个目标函数，使目标函数与输入的所有数据点最为接近，即利用最小二乘法所有数据点到曲线距离的平方和最小。假设已知有n个点 $\left(x_i,y_i\right)$ ( $i=1,2,\cdots ,n$ )，则拟合曲线的数学模型显示为：

$\begin{array}{cc}P\left(x\right)=a_0+a_{1}x+a_2{x}^2+\cdot \cdot \cdot +a_m{x}^m& \left(m<n\right)\end{array}$ (9)

使用最小二乘法求解参数，使得 $y_i$ 的偏差的平方和最小，即：

$R^2={\displaystyle \underset{i=1}{\overset{n}{\sum }}{\left[y_i-\left(a_0+a_{1}x+a_2{x}^2+\cdot \cdot \cdot +a_m{x}^m\right)\right]}^2}$ (10)

对等式右边求 $a_i$ 的偏导数：

$\{\begin{array}{l}\frac{\partial }{\partial a_0}{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\left[y_i-\left(a_0+a_{1}x+a_2{x}^2+\cdot \cdot \cdot +a_m{x}^m\right)\right]}^2=-2{\displaystyle \underset{i=1}{\overset{n}{\sum }}\left[y_i-\left(a_0+a_{1}x+a_2{x}^2+\cdot \cdot \cdot +a_m{x}^m\right)\right]=0}}\hfill \\ \frac{\partial }{\partial a_1}{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\left[y_i-\left(a_0+a_{1}x+a_2{x}^2+\cdot \cdot \cdot +a_m{x}^m\right)\right]}^2=-2{\displaystyle \underset{i=1}{\overset{n}{\sum }}\left[y_i-\left(a_0+a_{1}x+a_2{x}^2+\cdot \cdot \cdot +a_m{x}^m\right)\right]x=0}}\hfill \\ \frac{\partial }{\partial a_2}{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\left[y_i-\left(a_0+a_{1}x+a_2{x}^2+\cdot \cdot \cdot +a_m{x}^m\right)\right]}^2=-2{\displaystyle \underset{i=1}{\overset{n}{\sum }}\left[y_i-\left(a_0+a_{1}x+a_2{x}^2+\cdot \cdot \cdot +a_m{x}^m\right)\right]x^2=0}}\hfill \\ \cdots \hfill \\ \frac{\partial }{\partial a_m}{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\left[y_i-\left(a_0+a_{1}x+a_2{x}^2+\cdot \cdot \cdot +a_m{x}^m\right)\right]}^2=-2{\displaystyle \underset{i=1}{\overset{n}{\sum }}\left[y_i-\left(a_0+a_{1}x+a_2{x}^2+\cdot \cdot \cdot +a_m{x}^m\right)\right]x^m=0}}\hfill \end{array}$ (11)

将公式化简并整理为矩阵的形式：

$\left[\begin{array}{cccc}n& {\displaystyle \underset{i=1}{\overset{n}{\sum }}{x}_i}& \cdots & {\displaystyle \underset{i=1}{\overset{n}{\sum }}{x}_i^m}\\ {\displaystyle \underset{i=1}{\overset{n}{\sum }}{x}_i}& {\displaystyle \underset{i=1}{\overset{n}{\sum }}{x}_i{}^2}& \cdots & {\displaystyle \underset{i=1}{\overset{n}{\sum }}{x}_i{}^{m+1}}\\ ⋮& ⋮& \cdots & ⋮\\ {\displaystyle \underset{i=1}{\overset{n}{\sum }}{x}_i{}^m}& {\displaystyle \underset{i=1}{\overset{n}{\sum }}{x}_i{}^{m+1}}& \cdots & {\displaystyle \underset{i=1}{\overset{n}{\sum }}{x}_i{}^{2m}}\end{array}\right]\left[\begin{array}{l}{a}_0\\ a_1\\ ⋮\\ a_m\end{array}\right]=\left[\begin{array}{l}{\displaystyle \underset{i=1}{\overset{n}{\sum }}{y}_i}\\ {\displaystyle \underset{i=1}{\overset{n}{\sum }}{x}_i{y}_i}\\ ⋮\\ {\displaystyle \underset{i=1}{\overset{n}{\sum }}{x}_i{}^m{y}_i}\end{array}\right]$ (12)

上式可表示成 $X\ast A=Y$ ，则 $A=X^{-1}\ast Y$ ，如此便得到了系数矩阵A的值，从而得到拟合曲线的表达式。

5. 实验结果及分析

本文使用solidworks绘制带有孔洞的涡轮叶片，对数据进行3D打印，对打印后的模型进行扫描，如图3所示。

(a) (b)

Figure 3. Aeroengine turbine blade data. (a) Turbine Blade mechanical drawing; (b) 3D printing model

图3. 航空发动机叶片数据。(a) 叶片机械图；(b) 3D打印模型

通过对Alpha Shape算法的α值不断调整，实现了对涡轮叶片曲面点云的边缘提取。

Table 1. Different α bottom edge extraction data graph

表1. 不同α下边缘提取数据图

根据表1可以发现，当α = 0.6和α = 0.5时，提取的点数和提取时间有了较大幅度的增长，再结合两者提取的边缘图像来看，如图4所示。

图4(a)为α = 0.6时的边缘提取图，图4(b)为α = 0.5时的边缘提取图，两者对比不难发现，图4(b)的边缘提取由于α值过小导致提取失败，而图4(a)的边缘提取结果十分理想，所以采用α = 0.6作为边缘提取的参数。将通过Alpha Shape算法提取的平面边缘还原为三维曲面边缘，将三维曲面边缘使用欧几里德聚类将点云的内外边缘分割，从而获得孔洞区域的边缘点云，如图5所示。

(a) (b)

Figure 4. Plane edge extraction diagram. (a) Edge extraction image at α = 0.6; (b) edge extraction image at α = 0.5

图4. 平面边缘提取图。(a) α = 0.6时的边缘提取图；(b) α = 0.5时的边缘提取图

(a) (b)

Figure 5. Surface point cloud segmentation diagram. (a) Inner and outer edge segmentation diagram; (b) the edge of the segmented hole area

图5. 曲面点云分割图。(a) 内外边缘分割图；(b) 分割出来的孔洞区域边缘

通过观察线激光点云可以发现其特点：点云数据是由很多行线点云组成，而涡轮叶片的模型是同等曲率的曲面结构，所以可以通过完好曲线的线点云拟合一个曲线的多项式。由于我们的数据中的x，y轴是由扫描时的横向，纵向坐标组成，所以每条线点云的x，y坐标是一致的，所以只要将完整曲线的x，z轴的坐标代入用最小二乘法拟合多项式曲线，然后将孔洞区域的x轴坐标代入求解z轴坐标就可以修复孔洞区域的点云。

Table 2. Polynomial fitting analysis table

表2. 多项式拟合分析表

将提取的孔洞边缘按照行数切片。选取完整的点云线进行最小二乘多项式曲线拟合，从1次曲线到9次曲线依次拟合，对拟合后的多项式进行回归分析，根据拟合后的多项式的均方根误差以及拟合优度选取最优的多项式函数。根据表2可以看出，当多项式次数为8时，多项式曲线的拟合效果最好，得到多项式曲线参数为：

$\begin{array}{c}y=-9.763\times {10}^{-14}{x}^8-2.753\times {10}^{-11}{x}^7-3.31\times {10}^{-9}{x}^6-\text{2}\text{.093}\times {10}^{-7}{x}^5\\ -7.744\times {10}^{-6}{x}^4-2.128\times {10}^{-3}{x}^3-1.309\times {10}^{-2}{x}^2-0.85x-10.46\end{array}$

将孔洞区域的线点云的x，z轴坐标导入多项式曲线函数，获得修复后的曲线，如图6所示。

(a) (b)

Figure 6. Curve repair diagram. (a) Curve before repair; (b) curve after repair

图6. 曲线修复图。(a) 修复前的曲线；(b) 修复后的曲线

曲线拟合的结果如图6所示，图6(a)为修复前的曲线数据，图6(b)带有孔洞的线点云的x，z轴经过最小二乘法拟合后的曲线，通过图6(b)可以发现，橙色线条完美衔接了蓝色曲线的间断处。图7为涡轮叶片的叶背修复图，图中绿色部分点云为修复后的孔洞区域，通过修复后的可视化图像看，孔洞修复区域的点云和周围的点云贴合情况十分良好，因而孔洞修复的效果可以用作后续航空叶片的三维重建提供良好的基础。

Figure 7. Blade back repair drawing

图7. 叶背修复图

6. 结论

在叶片的使用过程中，叶片出现孔洞的情况在所难免，叶片的修复可极大降低维护成本。针对叶片的孔洞修复，本文使用基于线激光交叠式扫描的方式获取叶片的表面点云，使用统计滤波算法滤除扫描中产生的噪声，根据线激光点云由大量一行一行线点云组成的结构特点，提出了基于线激光的最小二乘曲线拟合的点云补全方法。补全方法步骤：首先使用基于Alpha Shape算法的曲面边缘提取方法提取航空叶片的内外边缘，再使用k近邻算法搜索的欧几里德点云聚类算法，通过欧几里德分割算法完成对使用Alpha Shape提取的叶片点云边缘的分割，最后根据获得的孔洞区域边缘得到每行缺陷曲线的数据，使用最小二乘拟合算法拟合叶片中完成的曲线，利用得到的多项式参数完成对缺陷曲线的补全，通过对一行行孔洞区域的曲线进行补全，完成叶片孔洞区域的修复。本方法是基于最小二乘法对线点云进行曲线拟合，由于最小二乘能够快速地拟合曲线，所以算法计算无需占用大量的算力和时间，是一种简单、快捷的点云孔洞修复方法。

参考文献