黄胜利 卜 昆 程云勇 周丽敏

西北工业大学现代设计与集成制造技术教育部重点实验室,西安,710072

0 引言

在涡轮叶片模具型腔反变形优化设计中,计算分析叶片测量数据相对于CAD模型的三维偏差至关重要。叶片测量数据的主要获取方法为三坐标测量和光学扫描测量。近年来,随着激光技术的发展,光学扫描仪获取数据的速度和精度都有显著的提高,应用越来越广泛,但其获得的海量数据给后续的数据处理、分析、计算带来很大的不便。由于测量数据是在测量设备坐标系下获得的,故在对测量数据进行偏差分析前,需把测量数据与CAD模型的坐标系统一起来,也就是实现测量数据与CAD模型的配准定位。

目前应用比较广泛的配准算法是Besl等[1]提出的迭代最近点(iterative closest point,ICP)算法。该方法通过最小二乘法构造配准的目标函数,具有收敛速度快、配准精度高的优点。基于ICP的改进算法[2]在含有自由曲面的工件配准中得到了广泛的应用。ICP算法要求待配准两模型的初始位置要足够接近,否则不能得到好的收敛结果,因此,使用ICP算法进行精确配准前,需要对两待配准模型进行预配准。武殿梁等[3]采用遗传算法实现了曲面之间的预配准,此方法比较耗时。吴锋等[4]基于轮廓的力矩主轴实现了医学图像的预配准,该方法的缺点是对数据的缺失比较敏感,需要整个物体的全部信息。刘晶[5]通过建立模型的包围盒实现三维测量模型与CAD模型之间的预配准,该方法对测量模型的完整性要求较高,且包围盒的建立有一定的难度。

光学测量的密集点云数据的配准计算非常耗时,为了缩短配准计算的时间,本文提出一种基于简化模型曲率计算的配准方法,首先对密集点云模型进行简化,估算简化模型的曲率,根据曲率提取对应特征面;然后寻找匹配点对从而计算出变换矩阵,把求出的变换矩阵应用于简化前的原始点云模型以实现模型的预配准;再通过奇异值分解和最近点迭代相结合(singular value decompo-sition and iterative closest point,SVD-ICP)算法[6-7]进行精确配准;最后以某型涡轮叶片的蜡模为例,验证方法的有效性和实用性。

1 点云模型简化

由于光学测量获取的数据量巨大,为了提高计算的速度需对点云模型进行适当的简化。图1所示分别为某型叶片及其蜡模的光学测量模型。目前,常用的点云简化方法有包围盒法、均匀网格法、聚类法、随机采样法、曲率采样法、均匀采样法等[8],它们各有优缺点和不同的适用范围。本文对点云模型进行简化的主要目的是提高计算的速度,因此采用均匀采样法。

图1 光学测量密集点云模型

2 曲率估算

曲率是曲面的重要几何特征,具有平移、旋转和缩放不变性,可以作为曲面特征识别的重要依据。本文采用邻域内拟合二次参数曲面的方法估算散乱点云的曲率[9]。

2.1 散乱数据参数化

对原始点云进行空间划分,计算点云中每个点pi的K 邻域,K邻域记为Knbhd(pi)。每个Knbhd(pi)可以拟合成一曲面且此曲面可以被一个局部基面逼近,该局部基面记为LBS(Knbhd(pi)),本文选择点pi处的切平面为局部基面。切平面可通过各点的邻近关系和最小二乘原理[10]求得。对于不同的模型,合理选取K的值(K一般取10～30),以保证邻域内凹凸一致,这样得到的切平面才能更好地逼近原始点集,使投影点在局部基面的参数化更好地反映点云的参数化。

把Knbhd(pi)中的各点在相应局部基面上的投影点记为PRO(Knbhd(pi)),投影点一般是随机分布的。为了参数化PRO(Knbhd(pi)),首先需要确定 u、v两个方向。在 PRO(Knbhd(pi))中 ,p′i是pi在局部基面上的投影点,求出距离p′i最远的点p′,将连接p′i、p′两点的直线方向作为u 方向,垂直于u方向的直线方向作为v方向。这里v方向可由u矢量和局部基面的法矢作向量积求得。然后将PRO(Knbhd(pi))中的每一点p′j(j=1,2,…,K)都与p′i相连,得到K 条有向线段,把这些有向线段分别与u矢量作点积,记为Dj,把得到的点积进行排序,最大值为Dmax,最小值为Dmin。由式(1)就可求得PRO(Knbhd(pi))中各点的u参数值:

其中,j=2,3,…,K+1且U1=0。同理可以求得各点的v参数值。把所求的PRO(Knbhd(pi))中各点的参数值作为Knbhd(pi)中对应点的参数,就实现了局部散乱数据的参数化。

2.2 拟合二次参数曲面

由2.1节中方法得到每个点pi的Knbhd(pi)的参数u、v,就可以建立二次参数曲面:

其中,A为3×3的系数矩阵。运用Yang[11]提出的二次参数曲面逼近法求得散乱数据点的二次参数曲面。首先引入新的矢量表达式:

则二次参数曲面可表示为

其中,系数a、b、c的下标为式(2)中A的行号和列号。采用最小二乘法使Knbhd(pi)上各点到二次参数曲面的欧几里德距离之和最小,可得

式中,X为邻域Knbhd(pi)的矩阵。

2.3 曲率计算

由于计算得到的曲面法矢方向一般是不协调的,即并不是所有的法矢都指向曲面的同一侧,因此有必要对求得的法矢进行归一化调整,具体方法参见文献[9]。

根据曲面的曲率性质可以得出高斯曲率:

平均曲率:

上述物理量可由曲面的第一、第二基本形式推出。

3 提取特征面

涡轮叶片主要由自由曲面和一些平面组成,而自由曲面的曲率特征没有一定的规律性,因此,曲面的提取有一定的难度。本文中,将测量模型和CAD模型上某一组对应平面作为寻找匹配点对的依据。高斯曲率I和平均曲率H均分别反映了曲面的凸凹形状,I和H的不同组合表现的是不同的几何形状,表1给出了几种常见的曲面形状与I、H组合的关系。

表1 点的局部曲面类型

由表1可以看出,当I、H 都为0时,曲面就是一平面,根据此条件可提取模型的某一平面为特征面,为寻找对应匹配点作好准备。

4 模型配准

4.1 预配准

为了分析铸件相对于其CAD设计模型的弯扭变形,首先需要实现测量模型与CAD模型的配准。运用SVD-ICP算法进行精确配准,要想获得较好的效果,首先需要实现测量模型与CAD模型的预配准。本文通过曲率提取特征面的方法实现预配准,具体方法如下:

(1)运用均匀采样法对原始测量点云{P}进行简化得到{P′}。

(2)计算简化后点云{P′}和CAD模型的重心(分别记为O1、O2),平移两模型,使两模型的重心与坐标原点O重合。

(3)运用第2节中所述方法估算简化后测量模型和CAD模型的高斯曲率K与平均曲率H。

(4)根据所求得的K和H提取两模型中相对应的特征面,分别记为L1、L2,并求取两特征面的

形心 P1、Q1。

(5)把P1、Q1沿所在特征面的法线方向移动一定距离 d(d值可适当选取,本文实验中 d=10mm),得点P2、Q2的计算公式:

式中,P1、P2、Q1、Q2为点 P1、P2、Q1、Q2的位置矢量;N1、N2分别为特征面的单位法矢。

(6)由三组对应点{O,P1,P2}和{O,Q1,Q2}求取旋转矩阵R和平移矩阵T,方法如下:以O点为坐标系原点,由单位矢量组 e1 、e2 、e3 与e′1 、e′2 、e′3按右手法则构建两个局部坐标系;设存在旋转变换矩阵R与平移变换矩阵T使两个局部坐标系重合,于是有

由式(13)可求得旋转变换矩阵:

待配准模型中的任一点Pi变换到对应点Qi的关系式为

把P1、Q1代入式(15)可得到平移变换公式:

(7)把求得的旋转矩阵和平移矩阵应用到简化前的原始测量点云模型上,得到变换后的测量模型,再按向量OO2进行整体平移就实现了测量模型与CAD模型的预配准。

4.2 精确配准

实现了测量模型与CAD模型的预配准后,就可以运用SVD-ICP算法进行精确配准。记测量点云数据为{pi},与其对应的CAD模型上的对应点集为{p′i}。配准目标函数表示为

对配准点集{pi}和{p′i}运用文献[7]中奇异值分解法求取变换矩阵R和T,并将变换矩阵作用到测量点云数据上,迭代地进行这个操作,直到使目标函数值满足设定的阈值为止。这样就实现了测量点云模型与CAD模型坐标系的统一,即实现了测量模型到CAD模型的精确配准。

5 实例验证

为了验证本文方法的有效性,以VC++6.0为平台,以VTK可视化类库为辅助工具编程实现了上述方法,并以某型涡轮叶片的蜡模为例进行了验证。试验中所用的测量点云模型是由ATOS流动式光学扫描仪获取的,其数据点数为102 908。试验所用的电脑配置为Intel(R)Core(TM)2 Duo CPU E6550@2.33GHz;内存为2.00GB。图2所示为蜡模测量模型与CAD模型的初始位置,图3所示为经上述方法计算得到的匹配点对。图4和图5分别为预配准和精配准的结果图。

图2 蜡模配准前初始位置

图3 对应匹配点对

图4 预配准结果

图5 精配准结果

运用本文方法对原始蜡模测量模型进行不同程度的简化,实现预配准和精确配准,并输出配准后的数据,通过Geomagic Studio 7.0对配准后的测量数据和CAD模型进行三维偏差计算,结果如图6所示。表2记录了不同简化程度下,配准计算所消耗的时间;表3记录了不同简化程度下对应的三维偏差结果。

图6 三维偏差结果

表2 不同简化程度下配准所耗时间对照表

表3 不同简化程度时三维偏差对照表

由表2可以看出,对原始测量点云数据进行不同程度的简化可以减少预配准所消耗的时间,精确配准时由于仍然是对原始测量点云数据进行计算,所以,此阶段所耗时间是相同的。由表3可知,与由原始点云计算的偏差相比,简化到不同点数时得到的三维偏差结果相差很小,以简化到12 864个点时的结果来看,与以原始点云作为控制点集计算得出的结果相比,其配准精度并没有受到影响,而其预配准阶段的时间从30s降到了4s,可见速度有了很大的提高。综合表2和表3的结果可知,对原始测量点云进行不同程度的简化,可以提高配准计算的速度,并且保证了配准的精度。

6 结语

本文提出一种基于简化模型的曲率计算配准方法。实例表明,该方法既保证了配准的精度又提高了配准的速度。但该方法只提高了预配准阶段的速度,并没有缩短精确配准阶段所耗的时间。

[1]Besl P J,MaKay N D.A Method of Registration of 3D Shapes[J].IEEE Transaction on Pattern Analysis and Machine Intelligence,1992,14(2):239-255.

[2]Li Y D,Gu P H.Free-form Surface Inspection Techniques State of the Artreview[J].Computer-Aided Design,2004,36(3):1395-1417.

[3]武殿梁,黄海量,丁玉成,等.基于遗传算法和最小二乘法的曲面匹配[J].航空学报,2002,23(3):285-288.

[4]吴锋,钱宗才,杭洽时,等.基于轮廓的力矩主轴法在医学图像配准中的应用[J].第四军医大学学报,2001,22(6):567-569.

[5]刘晶.叶片数字化检测中的模型配准技术及应用研究[D].西安:西北工业大学,2006.

[6]戴静兰,陈志杨,叶修梓.ICP算法在点云配准中的应用[J].中国图像图形学报,2007,12(3):517-521.

[7]张二虎,卞正中,张燕,等.基于ICP和 SVD的视网膜图像特征点配准算法[J].小型微型计算机系统,2004,25(10):1810-1813.

[8]戴静兰.海量点云预处理算法研究[D].杭州:浙江大学,2006.

[9]贺美芳,周来水,神会存.散乱点云数据的曲率估算及应用[J].南京航空航天大学学报,2005,37(4):515-519.

[10]张丽艳,周儒荣,蔡炜斌,等.海量测量数据简化技术研究[J].计算机辅助设计与图形学学报,2001,13(11):1019-1023.

[11]Yang M,Lee E.Segmentation of Measured Datausing a Parametric Quadric Surface Approximation[J].Computer-Aided Design,1999,31(7):449-457.