-
考虑到叶片形状的因素,将叶片大致分为叶盆、前缘、叶背和后缘等4个部分曲面,如图 1所示。鉴于每部分需要多次扫描获得数据,提出分阶段配准算法。先对这4个部分各自的多次数据进行自配准,即通过降维迭代最近点算法(dimensionality reduction iterative closest point,DRICP)来加速配准速度;然后每部分相邻之间采用改进代价函数的对称ICP算法(symmetric iterative closest point,SICP)进行配准,完成叶片全部的曲面配准。
-
通过研究,发现叶盆曲面和叶背曲面的各自数据集可以投影到某个平面后,与自身数据集不存在或者很少有重叠点。基于这个思路,提出一种快速降维配准算法,流程是先将各数据集依照某个平面投影降维,进行粗配准,然后计算点云的法线和曲率进行精配准。具体流程见下。
在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) 式中,(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]来采集叶片的深度信息,进而得到叶片点云数据。
航天发动机叶片具有高度扭曲性以及区域分段函数的特点。为了提高待配准点云的质量,提出动态规划扫描采集策略,保证在传感器有效测量范围内尽可能获取更多点云。因此,为了便于后续点云配准,将以实验平台的基准原点建立空间笛卡尔直角坐标系,扫描策略如图 4所示。图中红色框和绿色框代表前后两次扫描区域,箭头线段代表传感器运动轨迹。
考虑到叶片的扭曲程度以及传感器的有效量程,线结构光传感器仅沿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以内。
表 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算法在时间和误差有一定的均衡,使其优于前三者。
局部配准效果如图 7所示。可以看到,由于经典ICP算法对配准初始位置非常敏感,导致实验出现过匹配的现象,无法得到正确配准模型。GICP整体配准效果较好,但是边界处明显分层,有一定的夹角,有待进一步优化算法。NICP算法在主体部分配准较好,但是在顶部出现了明显的配准错误。自配准算法在重叠区域可以很好地完成配准任务,效果良好。
-
依照上面的配准实验后,可以得到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所示,可以方便地和原始模型进行差分运算,进行后续缺损部分的修复与测量。
面向航空损伤叶片点云的分阶段配准研究
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.
-
表 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 表 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 表 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 -