李景俊,芮执元,剡昌锋,王文斌,魏财斌
(1.兰州理工大学机电工程学院,兰州 730050;2.甘肃路桥飞宇交通设施有限责任公司,兰州 730050)
航空发动机叶片是航空飞行器的核心部件之一,其健康状况直接影响飞行器的安全性能和使用寿命[1]。涡轮叶片受载情况复杂,工作环境恶劣,容易造成疲劳损伤,而三维模型的逆向重建是修复受损叶片的关键技术之一[2]。缺陷孔洞的边界作为表达曲面的重要几何特征[3],边界特征点的快速、准确提取对于重建曲面模型具有重要意义[4]。
目前,航空发动机叶片缺陷孔洞的边界提取主要有两种方法,一种是在不确定数据结构的情况下基于三角网格提取缺陷孔洞边界[5],提取效果较好。但是对于海量点云,网格重构的复杂度高,占用内存空间大,耗时较长。另外,重构中不匹配的相接系数会产生错误的拓扑结构,导致重建模型结果失真或小孔洞识别失败,影响检测结果的准确性。另一种是基于三维点云的空间特征直接提取边界特征点。Nguyen等[6]采用区域网格增长算法确定目标网格边界,该算法效率较高,但是没有考虑噪声点对提取结果的影响。王小超等[7]利用协方差矩阵求得每一个采样点的法向度量,通过对阈值筛选后的特征三角形集进行子邻域划分并拟合每个子类点集的平面,判断采样点是否落在多个平面从而提取边界特征点。刘庆等[8]建立了点云的空间拓扑关系,基于K邻域内的采样点在微切平面上构建局部坐标系,以矢量合成的原理作为点云分布特性判别准则,实现边界特征点的提取。刘增艺等[9]采用K–D树建立了每个采样点及其K邻域点的空间拓扑关系,并基于K邻域点构造微切平面,通过判断微切平面上投影点之间的最大角度值是否大于角度阈值判别边界特征点。Qiu等[10]利用空间分割策略对点云模型的最小包围盒进行光栅化,建立了每个离散点的空间光栅索引和拓扑结构,在微切平面上根据邻域内投影点之间的夹角值筛选点云模型的边界特征点。王春香等[11]提出了一种针对密度不均、几何形状复杂点云模型的缺陷孔洞边界点识别方法,该方法采用相邻向量间角度序列的均值与比例因子的乘积作为最大角度阈值,对于非尖锐边界特征点的检测结果较好。Mineo等[12]提出了BPD算法判别边界特征点,该算法不需要设定任何阈值参数,通过局部邻域分辨率与邻域投影点在极坐标下的拓扑轨迹判别凹、凸型区域的边界特征点,适用性较好,但是运算复杂度较高。
对于几何形状不规则、采样点密度不均匀、曲率变化明显的航空发动机叶片等复杂形面的点云模型[13],使用单一判定准则难以完整地提取出多种类型的边界特征点,本文提出一种融合采样点的法向夹角判定准则和邻域场力矢量和判定准则判别边界特征点的方法,提高了边界特征点判别的准确度和适用性。为提高孔洞边界线的可连续性和光滑均匀程度,采用4次B样条拟合曲线对边界线进行平滑拟合,达到了较好的拟合效果。
本算法采用K–D树建立点云的空间拓扑结构以提高近邻点的搜寻效率,通过主成分分析法计算采样点的特征法向量,将采样点法向量与其邻域点法向量之间的夹角序列均值作为最大夹角阈值。基于K邻域点的场力矢量和准则,通过加权的方式获得采样点的特征值与判别阈值作为边界特征点的判定依据。为提高边界线的光滑均匀程度,采用4次B样条曲线拟合缺陷孔洞的边界线,算法流程图如图1所示。
数据采集卡收集的点云数据在空间上呈散乱分布,因此需要建立三维点云的空间拓扑结构。近邻域搜索方法主要有栅格划分法、八叉树搜索法、K–D树搜索法和R树搜索法等,其中K–D树搜索法具有搜索速度快,效率高等优点,因此采用K–D树对采样点数据进行空间结构划分,划分示意图如图2所示,可有效降低算法的复杂度,缩短算法的运行时间。
(1)法向夹角判定准则。
法向夹角判定准则根据相邻采样点的曲率特征变化实现对尖锐区域特征点的提取。本文基于K–D树算法,采用主成分分析法构造采样点pi及其邻域点qj的协方差矩阵E,通过求解该矩阵的特征值和特征向量得到采样点的特征法向量,
其中,E为协方差矩阵;K为邻域内的采样点数;qj为K邻域内的采样点;dj为K邻域的质心,。
对协方差矩阵E进行奇异值分解得到其特征值及特征向量,其中,最小特征值对应的特征向量Vmin即为采样点pi的特征法向量n0=(a,b,c)。
采样点与其邻域点之间的法向夹角可以反映采样点邻域曲面特征的变化情况。定义αi为任意采样点pi的法向量n0与其K邻域点qj法向量nj的夹角值,有
图1 算法流程图Fig.1 Flow chart of algorithm
在处理密度不均匀或边界中存在较多尖锐特征点的曲面时,可能存在边界特征点与较远的某一个非边界特征点形成较大夹角而造成特征点的误判。通过引入近邻域内两采样点之间的欧式距离作为调控系数以提高边界特征点判别的准确性和稳定性[14]。以采样点与其邻域点之间的法向夹角值代替常规两采样点之间的夹角值,新定义的局部邻域内采样点法向量之间的夹角值θi和最大角度阈值θt分别由式(3)和式(4)计算得出,
θi值的大小反映局部邻域曲面特征的尖锐程度,θi值越大表示该采样点的邻域曲面特征越尖锐,则该采样点属于边界特征点的概率越大。
(2)邻域场力矢量和判定准则。
邻域场力矢量和判定准则是基于邻域场的合力在同方向增强、反方向削弱的原理对平坦和低曲率区域进行边界特征点的判别。基于采样点pi的法向量方向,通过将K邻域内的任意一个采样点的坐标作为构建拟合曲面的约束条件,以K邻域内采样点作为局部平面参考依据构建拟合平面,该拟合平面的法向量为n0=(a,b,c),拟合平面和投影向量如图3所示,投影点处的坐标向量由式(5)计算得出,
图2 平面数据点的空间划分示意图Fig.2 Spatial division of plane data points
其中,v0=(x0,y0,z0)为平面上的任意采样点处的向量;vi=(xi,yi,zi)为平面外任意采样点处的向量;v′i=(x′i,y′i,z′i)为平面外任意采样点到拟合平面上的投影点处的向量。
对于拟合平面上投影点的分布,如果K邻域投影点在拟合平面上呈均匀分布,则判定该采样点为非边界特征点。如果其K邻域投影点在拟合平面上的分布偏向一侧,则判定该采样点为边界特征点[15]。如图4所示,定义采样点的投影点p′i为坐标原点,与距离p′i最远的投影点构成向量的方向作为X轴方向,以垂直X轴的方向为Y轴方向。通过将邻域采样点投影到拟合平面上得到邻域采样点的投影点坐标为q′j=(q′jX,q′jY)。
图3 采样点处拟合平面和投影向量示意图Fig.3 Fitting tangent plane and projection vector
如图4所示,采样点的投影点所构成的向量在拟合平面上一般呈无序分布,为计算简便,根据式(6)将投影向量进行单位化处理,单位化后的投影向量如图5所示,令单位化后K邻域内采样点的投影点坐标为q″j(Xj,Yj)。
其中,p′j q″j为单位化后的采样点与其K邻域内采样点之间的投影向量;p′i为采样点的投影点,q″j(j=0,1,2,…,K–1且j≠i)为单位化后的K邻域内采样点的投影点。
赋予向量p′j q″j一个作用力Fij,力的方向与其向量方向相同,力的大小等于|p′j q″j|,通过计算K邻域内投影向量的合力在X轴和Y轴上的分力和来判别边界特征点,分力分别由式(7)和(8)计算得到,
图4 K邻域内采样点的投影示意图Fig.4 Projection of sampling points in K neighborhood
图5 单位化后的投影向量Fig.5 Projection vector after normalization
邻域投影向量的合力在X轴或Y轴方向上的分力值越大,则该方向上的采样点属于边界特征点的概率越大。因此取邻域投影向量的合力在X轴或Y轴方向上较大分力值作为邻域场力和,记作Fi。
为提高缺陷孔洞边界特征点提取的精准度,本算法融合上述两判定准则的优点,引入μ1和μ2分别作为法向夹角判定准则和邻域场力矢量和判定准则的调和系数,采用加权方式实现边界特征点的提取,可获得更加全面、准确的边界特征点提取结果。根据求得的θi和Fi,分别由式(9)和式(10)计算得出数据采样点的特征值ti和判别阈值T:
其中,μ1+μ2=1,F为邻域内每个采样点在同方向的场力和。
其中,η是特征点数控制系数;N为初始采样点数。
当ti≥T时,则判定该采样点为边界特征点。为提高算法的自适应性和准确性,调和系数μ1和μ2的比值关系根据式(11)计算得出:
其中,K1为K邻域内θi≥θt的个数;K2为K邻域内|Fi|≥|Fe|的个数;Fe为采样点的投影点与其邻域内所有投影点作用力的平均值,
对于大多数复杂孔洞工件,其边界多为复杂曲线,常规的3次B样条拟合的边界线存在较多的锯齿状边缘,不能较好地满足复杂几何形状边界线拟合的需求[16]。本算法基于欧式距离将相邻的边界特征点建立一个元胞数组,根据由每个元胞数组构成的采样点集,采用4次B样条拟合曲线可有效地提高算法的精准度和曲线的光滑度。
4次B样条曲线至少需要5个控制点,因此在采样点数据的首尾两端分别添加两个补充向量,为了使数据结构更具有周期性,将原始数据中末端和首端两个采样点分别作为数据首尾两端的补充向量,补充后的数据向量为,
具体计算步骤如下。
Step(1):将补充向量作为拟合曲线的初始控制点,第k次迭代的控制点记为p<k>,令p<k>=p*。
Step(2):采用4次B样条矩阵形式,则第k次迭代后的近似B样条曲线s(t)
其中,s(t)
Step(3):记p
Step(4):根据迭代误差公式e=p*–p*
Step(5):定义新的控制点p
试验采用20MHz超声换能器,六自由度的机械手(STAUBLI TX90L)夹持叶片以垂直叶片法向量的位姿等间距地移动,叶片的扫查装置如图6所示。通过数据采集卡(采集频率:250MHz)实现航空叶片的原始数据采集,原始数据经过闸门设置和信号预处理后获得超声C扫点云数据。
图6 叶片的超声扫查装置实物图Fig.6 Ultrasonic scanning device of blade
叶片A是一块预制多处缺陷孔洞的航空发动机叶片,如图7(a)所示,缺陷孔洞的几何形状不规则、面积大小不同且位置分布散乱。叶片B为一块预制多种尺寸平底孔的航空发动机叶片,实物如图7(b)所示。表1为两块叶片点云数据的详细信息,表2为叶片B预制平底孔的尺寸信息。
边界特征点提取试验所使用的计算机的配置参数为:Intel(R) Core i7–7700 HQ,CPU 2.8GHZ,内存8GB,Windows 10专业版64位操作系统。
基于上述的理论基础,在具有多个复杂形状缺陷孔洞的航空发动机叶片点云模型上进行缺陷孔洞边界特征点的提取试验。分别使用单一判定准则[9]判别方法和本文所提出的算法进行试验对比分析,试验结果对比如表3所示。
图7 叶片实物图Fig.7 Blade photos
研究表明,邻域点数K的选取对边界特征点判别结果影响较大,当K值较小时不能准确反映该邻域点的分布情况,提取结果中存在较多的误判点,缺陷孔洞边界不清晰;K值选取较大则会覆盖小孔洞,导致小孔洞识别失败。K值的选取需要根据采样点的密度大小和噪声点的数量确定,本试验的两组数据选用K=20均获得较为满意的提取效果。另外,特征点数控制系数η用于调节判定阈值的大小进而控制边界特征点的数量,本试验选用η=0.7获得较好的提取结果。
通过对比表3试验结果和图8试验效果,使用单一判定准则难以完整地提取出复杂缺陷孔洞的边界特征点,而本算法可有效地提取更多的边界特征点,提高了边界特征点提取的准确度。对比图9试验效果,常规的边界点判别方法针对0.2mm较小孔洞边界特征点提取效果不佳,而本文采用的加权判断准则可以较完整地提取出小孔洞的边界特征点。引入两采样点之间的欧式距离作为调控系数,降低了远距离点对特征点判别的影响,对于密度分布均匀或非均匀的点云模型边界特征点提取均具有一定的适用性。对比由3次B样条拟合的边界线效果与本文采用的4次B样条拟合的边界线效果(图10和11),发现3次B样条拟合的边界线中产生较多锯齿状边缘,边界线不够平滑,而本文提出的曲线拟合算法达到较好的平滑效果,提升了复杂形面边界线的可连续性和光滑均匀程度,对于缺陷叶片的模型重建具有一定的实际应用价值。
(1)对于航空发动机叶片等复杂几何形面,提出一种融合法向夹角判定准则和邻域场力矢量和判定准则提取曲面模型中边界特征点的方法,试验结果表明,该方法可有效提取出多种复杂孔洞模型的边界特征点。
表1 超声C扫描点云数据参数Table 1 Specifications of ultrasonic C-scan point cloud
表2 叶片B预制缺陷孔洞的尺寸Table 2 Pre-fabricated defect specifications of blade B
表3 试验结果对比Table 3 Comparison of experimental results
图8 叶片A孔洞边界点提取效果图Fig.8 Boundary point cloud extraction of blade A
图9 叶片B孔洞边界点提取效果图Fig.9 Boundary point cloud extraction of blade B
图10 叶片A缺陷孔洞边界点拟合效果对比Fig.10 Comparison of fitting results of boundary points of defect hole in blade A
图11 叶片B缺陷孔洞边界点拟合效果对比Fig.11 Comparison of fitting results of boundary points of defect hole in blade B
(2)针对采样点与邻域点之间夹角过大而造成特征点误判的问题,通过引入两采样点之间的欧式距离作为调控系数,以采样点法向量与其邻域点法向量的夹角代替常规的采样点与其邻域点的夹角值,提高了边界特征点判别的准确性和适用性。