网站地图
-
随着我国航空事业的快速发展，航空发动机的研究愈发重要，而叶片作为高性能航空发动机的核心工作部件，极为重要。叶片在高温、高压、高摩擦的环境中高速旋转，承受较大负荷，损伤率高，且其造价昂贵，一般多采用修复技术以延长使用寿命。叶片的工作环境要求叶片在修复过程中必须采用高精度的检测方法以保证叶片修复可以达标[1]。因此，建立高精度的损伤叶片模型具有重要的实际意义。随着测量技术的发展，采用激光三角法原理[2]的非接触检测方法开始在航空叶片的损伤检测领域中广泛使用。而其数据处理算法还有待优化，许多学者对此展开了研究。
ZHOU等人[3]针对直纹叶片提出了线性特征构造，通过对边界点云数据进行非刚性配准和等分数据来快速得到叶片曲面，该方法具有一定局限性。MA等人[4]针对风电叶片特征点少的问题提出了基于直线特征的配准方法，提高了叶片重建的速度，但是难以得到叶片3维模型。SONG等人[5]提出采用移动最小二乘算法对叶片点云数据曲面重建，由于叶片自身是自由曲面，该方法获得的曲面较为平滑，但是处理过程可能形成奇异方程造成错误配准。ZHANG等人[6]通过约束最小二乘的方法对叶片边缘数据进行拟合建模，采用了以前缘为中心逐次拟合的方法增强算法的鲁棒性，可以提高建模精度，但是需要提前已知叶片模型边缘形状。
BESL等人[7]提出迭代最近点算法(iterative closest point algorithm，ICP)，该算法对点云的位置进行变化以达到配准的目的，但是对初始位置依赖严重，且容易陷入局部最优。LIU等人[8]提出将ICP算法用在叶型点云数据的配准中，使用四元素法进行求解，有一定的效果。LIN等人[9]提出3次样条的优化搜索策略来改进ICP算法对叶片的截面曲线数据进行配准，减少了迭代次数，但还是严重依赖初始配准位置。SEGAL等人[10]采用最近分布到分布的对应模型，在经典ICP算法里将代价函数中附加了一个高斯概率模型，增加了计算的复杂度。SERAFIN等人[11]提出基于法线特征与曲率的改进迭代最近点配准算法，对初始位置有一定的改进。CHEN等人[12]采用天牛须改进的粒子群算法以分布熵为寻优标准的方法来加速点云的配准，该方法还有待进一步优化。LIU等人[13]在参考文献[14]的基础上提出了混合稀疏的ICP算法，采用混合正则项作为配准函数，结合交替乘子法进行配准，提高了配准效率，但是对低重叠率叶片数据的配准有待优化。PROKOP等人[15]提出基于线特征的低重叠率点云配准，该方法非常依赖预处理和线特征提取的结果。
针对配准初始位置的优化，DAI等人[16]提出先采用主方向贴合的方法进行点云的粗配准，通过计算两块点云重心进行简单平移，该方法对形状差异较大的点云提供的初始位置并不理想。SHI等人[17]采用主成分分析(principal component analysis，PCA)的改进ICP策略，与经典ICP算法相比，明显缩短了整体配准时间，但配准时间还是很长。ZHAO等人[18]通过采样一致性进行初始位置的变换，进一步提高了后续配准的精度。LU等人[19]通过提取点云尺度不变特征来增加局部特征的权重，以此使用超四点快速鲁棒匹配算法[20]进行粗配准，可以提供较好的初始位置。
本文作者基于以上对叶片数据处理算法的研究，结合叶片自身特点，提出一种分阶段配准方法，将叶片分为4个部分曲面，对每部分曲面进行自配准，并对不同部分曲面进行完全配准。自配准先使用降维的方法进行粗配准，再使用改进的迭代最近点算法进行精配准；完全配准在前面的基础上再次配准，获得叶片完整曲面。搭建了基于线激光的四维度扫描系统实验平台，实验证明在满足精度的前提下，本文中的算法提高了叶片检测效率，得到较高精度的叶片3维模型。
-
考虑到叶片形状的因素，将叶片大致分为叶盆、前缘、叶背和后缘等4个部分曲面，如图 1所示。鉴于每部分需要多次扫描获得数据，提出分阶段配准算法。先对这4个部分各自的多次数据进行自配准，即通过降维迭代最近点算法(dimensionality reduction iterative closest point，DRICP)来加速配准速度；然后每部分相邻之间采用改进代价函数的对称ICP算法(symmetric iterative closest point，SICP)进行配准，完成叶片全部的曲面配准。
图 1 叶片曲面分布
Figure 1. Distribution of the blade surface
-
通过研究，发现叶盆曲面和叶背曲面的各自数据集可以投影到某个平面后，与自身数据集不存在或者很少有重叠点。基于这个思路，提出一种快速降维配准算法，流程是先将各数据集依照某个平面投影降维，进行粗配准，然后计算点云的法线和曲率进行精配准。具体流程见下。
在3维坐标系内选取一个平面M，使源点云集和目标点云集P和Q投影到该平面，得到新的自配准点云集Pz和Qz，且平面内无重叠点。为简便起见，取3维坐标系的原点投影到平面M的点O′为平面2维坐标系的原点，同时取另外一个轴线基向量做平面M的投影，可以此建立完整2维坐标系，如图 2所示，得到单位投影算子$\nabla$:
$ \nabla=\left[\begin{array}{lll} x_{\mathrm{h}} & y_{\mathrm{h}} & z_{\mathrm{h}} \\ x_1 & y_1 & z_1 \\ x_{\mathrm{w}} & y_{\mathrm{w}} & z_{\mathrm{w}} \end{array}\right] $
(1) 
图 2 坐标系转换
Figure 2. Coordinate system conversion
式中，(xn yn zn)(n=h，l，w)为投影算子在3维坐标下的对应值。
将原坐标系下的某点D(xD，yD，zD)投影变换到平面后的点d(xd，yd)，坐标转换公式如下：
$ \left(h_d, x_d, y_d\right)=\nabla \cdot\left[\begin{array}{c} x_D \\ y_D \\ z_D \end{array}\right] $
(2) 式中，hd为点D(xD，yD，zD)到平面M的距离。
计算自配准点云集Pz和Qz去质心化后的点云集Po和Qo，计算表达式为：
$ \left\{\begin{array}{l} p_{\mathrm{o}, i}=p_{\mathrm{z}, i}-\frac{1}{N_{\boldsymbol{P}_{\mathrm{z}}}} \sum\limits_{i=1}^{N_{\boldsymbol{P}_{\mathrm{z}}}} p_{\mathrm{z}, i} \\ q_{\mathrm{o}, i}=q_{\mathrm{z}, i}-\frac{1}{N_{\boldsymbol{Q}_{\mathrm{z}}}} \sum\limits_{i=1}^{N_{Q_{\mathrm{z}}}} q_{\mathrm{z}, i} \end{array}\right. $
(3) 式中，NPz，NQz为点云集Pz，Qz的点云数量，pz，i，qz，i分别为点云集Pz和Qz的点云，po，i，qo，i分别为去质心化的点云集Po和Qo中与之对应的点云。
寻找两个点云集Po，Qo的点对应关系。当点云集Po中某点po，i=(xpo，i，ypo，i)与点云集Qo中某点qo，i=(xqo，i，yqo，i)对应，取二者欧氏距离为评判标准，则对应点计算公式如下：
$ \begin{gathered} \left(p_{0, i}, q_{0, i}\right)= \\ \arg \min _i\left(\sqrt{\left(x_{q_{0, i}}-x_{p_{0, i}}\right)^2+\left(y_{q_{0, i}}-y_{p_{0, i}}\right)^2}\right) \end{gathered} $
(4) 计算平均欧氏距离，作为配准误差标准：
$ E_{\mathrm{o}}\left(\boldsymbol{R}_{\mathrm{o}}, \boldsymbol{T}_{\mathrm{o}}\right)=\frac{1}{N_{\mathrm{o}}} \sum\limits_{i=1}^{N_{\mathrm{o}}}\left\|q_{\mathrm{o}, i}-\left(\boldsymbol{R}_{\mathrm{o}} p_{\mathrm{o}, i}+\boldsymbol{T}_{\mathrm{o}}\right)\right\|_2 $
(5) 式中，No为点云集Po和Qo中具有对应关系的点对数量，为点云集Po的旋转矩阵，To为点云集Po的平移矩阵，θo为点云集Po旋转角度。通过对其展开化简可得:
$ \begin{gathered} E_{\mathrm{o}}\left(\boldsymbol{R}_{\mathrm{o}}, \boldsymbol{T}_{\mathrm{o}}\right)=\frac{1}{N_{\mathrm{o}}}\left(\sum\limits_{i=1}^{N_{\mathrm{o}}}\left|q_{\mathrm{o}, i}\right|^2+\right. \\ \left.\sum\limits_{i=1}^{N_{\mathrm{o}}}\left|\boldsymbol{R}_{\mathrm{o}} p_{\mathrm{o}, i}\right|^2-2 \sum\limits_{i=1}^{N_{\mathrm{o}}}\left|q_{\mathrm{o}, i} \boldsymbol{R}_{\mathrm{o}} p_{\mathrm{o}, i}\right|\right) \end{gathered} $
(6) 对其求导展开后可得：
$ \left\{\begin{array}{l} \theta_{\mathrm{o}}=\arctan \left[\frac{\sum\limits_{i=1}^{N_\mathrm{o}}\left(x_{q_{0, i}} y_{p_{0, i}}-y_{q_{\mathrm{o}, i}} x_{p_{0, i}}\right)}{\sum\limits_{i=1}^{N_\mathrm{o}}\left(x_{q_\mathrm{o}, i} x_{p_{0, i}}+y_{q_{\mathrm{o}, i}} y_{p_{0, i}}\right)}\right] \\ \boldsymbol{T}_{\mathrm{o}}=\frac{1}{N_{\mathrm{o}}} \sum\limits_{i=1}^{N_{\mathrm{o}}} q_{\mathrm{z}, i}-\boldsymbol{R}_{\mathrm{o}} \frac{1}{N_{\mathrm{o}}} \sum\limits_{i=1}^{N_{\mathrm{o}}} p_{\mathrm{z}, i} \end{array}\right. $
(7) 由(7)式求得降维后的旋转矩阵和平移矩阵，和点云集Pz一起代入(5)式计算误差。当误差小于设定阈值时，用原坐标系的旋转矩阵和平移矩阵，对点云集Pz进行反向变换，得到最新的自配准的点云集Pc。在这个阶段，点云集Q不做变换。计算点云集Q，Pc的各点法线特征和曲率特征，根据点云集的特征，采用图索引的方式，搜索二者正确的对应点云集合。
由对应点云集计算出旋转矩阵Rc和平移矩阵Tc，通过变换矩阵计算点云集Q，Pc的配准误差Ec(Rc，Tc)：
$ \begin{gathered} E_{\mathrm{c}}\left(\boldsymbol{R}_{\mathrm{c}}, \boldsymbol{T}_{\mathrm{c}}\right)= \\ \frac{1}{N_{\mathrm{c}}}\left[\sum\limits_i^{N_{\mathrm{c}}} \omega_{\mathrm{c}, i}\left\|\boldsymbol{\eta}_{q_i}\left(\boldsymbol{R}_{\mathrm{c}} q_i+\boldsymbol{T}_{\mathrm{c}}-p_{\mathrm{c}, i}\right)\right\|^2\right] \end{gathered} $
(8) 式中，$\boldsymbol{P}_{\mathrm{c}}=\left\{p_{\mathrm{c}, i} \mid i=1, 2, \cdots, N_{\mathrm{c}}\right\}, \boldsymbol{Q}=\left\{q_i \mid i=1, 2, \cdots, N_{\mathrm{c}}\right\}, $ ηqi为点qi邻域的表面法线。
$ \omega_{\mathrm{c}, i}=\left\{\begin{array}{l} 1, \left(\left\|p_{\mathrm{c}, i}-\boldsymbol{R}_{\mathrm{c}} q_i-\boldsymbol{T}_{\mathrm{c}}\right\| \leqslant d_{\max }\right) \\ 0, \left(\left\|p_{\mathrm{c}, i}-\boldsymbol{R}_{\mathrm{c}} q_i-\boldsymbol{T}_{\mathrm{c}}\right\|>d_{\max }\right) \end{array}\right. $
(9) 式中，dmax为设定的最大匹配距离阈值。
当误差小于设定值或者超出迭代次数时，停止迭代，完成自配准。
-
考虑到完成自配准后的点云数据很难再通过投影到某一平面来快速完成粗配准，为了将不同视角的点云整合到同一空间坐标系，将优化ICP算法的代价函数进行最终配准。
迭代最近点算法是目前应用最广泛的精确拼接算法，算法在每次迭代中缩小源点云到目标点云中最近点的欧氏距离，通过不断迭代完成配准。ICP算法操作简单，但是对配准初始位姿要求较高。通过改进ICP算法对不同曲面的配准，使其在初始位置不完美的情况下，采用点云的曲率法线快速完成配准，步骤如下：(1)对完全配准的点云集Pw，Qw进行对应点匹配; (2)计算点云的法向量npw，i，nqw，i; (3)计算对应点的3维旋转矩阵和平移矩阵Rw，Tw，并对源点云进行位置变换; (4)根据以下标准计算完全配准的配准误差Ew：
$ \begin{gathered} E_{\mathrm{w}}=\frac{1}{N_{\mathrm{w}}} \sum\limits_i^{N_{\mathrm{w}}}\left[\left(\boldsymbol{R}_{\mathrm{w}} p_{\mathrm{w}, i}-\boldsymbol{R}_{\mathrm{w}}{ }^{-1} q_{\mathrm{w}, i}+\boldsymbol{T}_{\mathrm{w}}\right) \times\right. \\ \left.\left(\boldsymbol{\eta}_{p_{\mathrm{w}, i}}+\boldsymbol{\eta}_{q_{\mathrm{w}, i}}\right)\right]^2 \end{gathered} $
(10) 式中，$\boldsymbol{P}_{\mathrm{w}}=\left\{p_{\mathrm{w}, i} \mid i=1, 2, \cdots, N_{\mathrm{w}}\right\}, \boldsymbol{Q}_{\mathrm{w}}=\left\{q_{\mathrm{w}, i} \mid i=1, 2, \cdots, \right.$Nw}，Nw为点云集Pw，Qw中认确匹配的对应点数量$\boldsymbol{\eta}_{p_{\mathrm{w}, i}}, \boldsymbol{\eta}_{q_{\mathrm{w}, i}} \text { 为 } p_{\mathrm{w}, i}, q_{\mathrm{w}, i}$各自的邻域表面的法线; (5)当误差小于设定值或达到迭代设置次数时停止计算并输出结果，否则更新点云Pw，Qw位置坐标数据，返回到步骤(1)继续计算。
-
搭建3维重建系统实验平台如图 3所示。系统利用三角测距法[2]来采集叶片的深度信息，进而得到叶片点云数据。
图 3 系统实验平台
Figure 3. System experimental platform
航天发动机叶片具有高度扭曲性以及区域分段函数的特点。为了提高待配准点云的质量，提出动态规划扫描采集策略，保证在传感器有效测量范围内尽可能获取更多点云。因此，为了便于后续点云配准，将以实验平台的基准原点建立空间笛卡尔直角坐标系，扫描策略如图 4所示。图中红色框和绿色框代表前后两次扫描区域，箭头线段代表传感器运动轨迹。
图 4 叶片扫描策略
Figure 4. Blade scanning strategy
考虑到叶片的扭曲程度以及传感器的有效量程，线结构光传感器仅沿yOz平面运动无法获取完整的叶片数据，存在许多无效值。为此提出在沿图 4的规划路径进行数据采集的同时，需要根据采集数据的质量进行动态改变x轴坐标。
同时，当相邻两次的采样区间可匹配特征较少时，会导致后续配准无法完成。为此提出了采样区间的重叠比δ，公式如下：
$ \delta=\frac{m-y_{i+1}+y_i}{m} $
(11) 式中，yi+1和yi为传感器前后两次移动的y轴坐标值，m为传感器的有效量程。
若当重叠比δ取值较小，则前后两次数据重叠特征过少，容易配准失败，当重叠比δ取值较大时，会增加后续配准的内存消耗和CPU处理时间，因此选取了重叠比为0.5来平衡配准速度和配准精度。
-
为获取系统测量误差，设计了一组0级的标准量块测量实验，如图 5所示。该标准量块宽度为8 mm，偏差为0.05 μm，可将宽度视为真值。控制传感器进行采集量块数据，记录量块宽度左右两侧位置坐标数据，如表 1所示。由表 1可知，系统测量误差在0.010 mm以内。
图 5 量块测量实验
Figure 5. Gauge block measurement experiment
表 1 测量标准量块宽度(单位：mm)
Table 1. Measuring standard gauge block width/mm
1 2 3 4 5 left edge 10.835 10.869 10.771 10.768 10.735 right edge 18.843 18.876 18.777 18.777 18.744 gauge block width 8.008 8.007 8.006 8.009 8.009 deviation 0.008 0.007 0.006 0.009 0.009 -
以叶盆曲面为例，选取重叠比为0.5进行数据采集，通过对叶盆的6次数据分别进行两两配准，本文中的DRICP算法与ICP算法、广义迭代最近点(generalized iterative nearest point algorithm，GICP)算法、法线迭代最近点(normal iterative nearest point algorithm，NICP)等参考算法进行对比实验，结果如表 2所示。由表 2可知，在精度近似相同的情况下，自配准算法所需时间在配准叶盆曲面的6次实验中均优于对比算法。
表 2 叶盆曲面配准实验时间(单位：s)
Table 2. Leaf basin surface registration experiment time/s
registration algorithm 1 2 3 4 5 6 ICP 5.626 9.066 7.178 8.757 10.598 5.949 GICP 48.892 106.937 67.453 102.838 134.986 68.009 NICP 2.021 6.656 2.239 9.871 2.978 1.713 DRICP 0.532 0.512 0.823 0.705 0.970 0.970 配准实验误差如图 6所示。经典ICP算法配准的出现错误匹配，导致误差过小；而GICP和NICP在第1次和最后1次的配准中波动较大，表明二者在曲面特征很少时配准过程还有待优化；GICP算法在其余4次实验中误差均小于本文中的算法，但是结合表 1来看，所耗时间过长；本文中的DRICP算法在时间和误差有一定的均衡，使其优于前三者。
图 6 配准误差曲线
Figure 6. Registration error curve
局部配准效果如图 7所示。可以看到，由于经典ICP算法对配准初始位置非常敏感，导致实验出现过匹配的现象，无法得到正确配准模型。GICP整体配准效果较好，但是边界处明显分层，有一定的夹角，有待进一步优化算法。NICP算法在主体部分配准较好，但是在顶部出现了明显的配准错误。自配准算法在重叠区域可以很好地完成配准任务，效果良好。
图 7 叶盆局部数据配准曲面效果
Figure 7. Effect of the local registration surface of the leaf basin
-
依照上面的配准实验后，可以得到4个叶片扫描曲面。首先通过PCA粗配准，得到一个较为良好的初始位置，再使用改进ICP算法进行配准，实验结果如表 3所示。从表 3可见，优化ICP算法的误差函数可以提高算法的准确性。
表 3 完全配准实验结果
Table 3. Experimental results of full registration
algorithm test 1 2 3 GICP time/s 5.223 4.741 5.661 error/μm 1.3 28.0 1.6 NICP time/s 4.885 4.613 3.602 error/μm 1.9 61.3 3.1 SICP time/s 4.69 4.224 5.374 error/μm 1.4 26.0 1.5 依次通过3次配准可以得到叶片完整模型，如图 8a所示，配准效果良好。导入到专业软件可得到叶片重建模型，如图 8b所示，可以方便地和原始模型进行差分运算，进行后续缺损部分的修复与测量。
图 8 处理后的航空发动机损伤叶片
Figure 8. Damaged blade of aero-engine after treatment
-
针对以非接触式测量方式采集的叶片数据处理算法不同导致模型精度偏差的问题，提出一种分阶段配准算法来优化叶片点云数据的处理，通过搭建线激光叶片3维扫描系统，对待检测航空叶片进行3维扫描，得到其完整的点云数据，继而采用本文中的算法配准重建得到叶片的3维模型，与现有算法相比，在精度相近的情况下，提高了测量速度，对后续叶片的增材修复以及产品的升级，具有重要意义和前景价值。
面向航空损伤叶片点云的分阶段配准研究
Research on staged registration of point clouds for aeronautical damaged blades
-
摘要: 为了提高非接触式测量的数据处理精度, 采用一种分阶段配准的方法, 先将缺损叶片分为4个部分, 采用自配准算法对每部分进行配准; 再对相邻两部分采用改进的完全配准算法进行整体配准。结果表明, 自配准算法与传统算法相比, 在配准误差均小于0.005 mm的前提下, 配准时间可以缩短到1 s以内; 完全配准与传统算法相比, 速度较快, 并通过0级的标准量块测量实验得出系统的测量误差小于0.010 mm, 满足叶片测量的精度要求。该分阶段配准方法对测量航空叶片具有一定的应用价值。Abstract: In order to improve the accuracy of data processing, a staged registration method was proposed. First, the defective blades were divided into four parts, and the self-registration algorithm was used to register each part; then the improved complete registration algorithm was used for the global registration of the adjacent two parts. The experimental results show that compared with the traditional algorithm, the registration time can be shorten to less than 1 s under the premise that the registration error is less than 0.005 mm by using the self-registration algorithm; the complete registration is faster than the traditional algorithm. The measurement error of the system is less than 0.010 mm by measuring the 0-level standard gauge block, which meets the accuracy precision requirements of blade measurement. The staged registration method has certain application value for the measurement of aviation blades.
-
图 8 处理后的航空发动机损伤叶片
a—配准后的叶片 b—重建后的叶片
Figure 8. Damaged blade of aero-engine after treatment
a—registered blade b—reconstructed blade
表 1 测量标准量块宽度(单位：mm)
Table 1. Measuring standard gauge block width/mm
1 2 3 4 5 left edge 10.835 10.869 10.771 10.768 10.735 right edge 18.843 18.876 18.777 18.777 18.744 gauge block width 8.008 8.007 8.006 8.009 8.009 deviation 0.008 0.007 0.006 0.009 0.009 
下载: 导出CSV
表 2 叶盆曲面配准实验时间(单位：s)
Table 2. Leaf basin surface registration experiment time/s
registration algorithm 1 2 3 4 5 6 ICP 5.626 9.066 7.178 8.757 10.598 5.949 GICP 48.892 106.937 67.453 102.838 134.986 68.009 NICP 2.021 6.656 2.239 9.871 2.978 1.713 DRICP 0.532 0.512 0.823 0.705 0.970 0.970
下载: 导出CSV
表 3 完全配准实验结果
Table 3. Experimental results of full registration
algorithm test 1 2 3 GICP time/s 5.223 4.741 5.661 error/μm 1.3 28.0 1.6 NICP time/s 4.885 4.613 3.602 error/μm 1.9 61.3 3.1 SICP time/s 4.69 4.224 5.374 error/μm 1.4 26.0 1.5
下载: 导出CSV
-
-
点击查看大图
图(8) / 表(3)
计量
- 文章访问数: 1052
- HTML全文浏览量: 571
- PDF下载量: 8
- 被引次数: 0
