CN110717276B - 基于工业ct扫描的异型气膜孔几何结构检测与评定方法 - Google Patents

Abstract

本发明提供了一种基于工业CT扫描的异型气膜孔几何特征检测与评定方法。利用工业CT扫描得到异型孔轮廓三维点云坐标，通过剔除非轮廓点、坐标变化、截面积分、三维配准等运算处理，分别提取获得异型孔几何尺寸、形状与轮廓度等评定参数。本发明方法测量精度高、实现简单，弥补了异型孔检测与评定标准的缺失，能够提高异型孔几何精度检测与评定的效率与精度。

Description

基于工业CT扫描的异型气膜孔几何结构检测与评定方法

技术领域

本发明属检测技术领域，具体涉及一种基于工业CT扫描的异型气膜孔几何结构检测与评定方法。

背景技术

发明专利：CN109974648A公开了一种基于坐标值的微小孔几何精度评定方法。该方法基于三维表面轮廓仪测量得到圆孔出入口平面的几何坐标，提取孔周轮廓点坐标，利用最小二乘法，最小区域法计算微小孔直径、圆度误差。该方法是针对圆柱型孔，利用圆柱孔出入口直径、圆度和锥度等表征参数对其几何结构进行评定。然而对于形状较为复杂的三维结构如异型孔，仅利用出入口平面的尺寸和轮廓度误差等表征参数，难以准确的反映异型孔的几何结构特征。

气膜冷却孔主要应用于航空发动机涡轮叶片中，气膜孔是在涡轮叶片前缘，叶身等位置加工直径约为0.2-0.5mm小孔，从榫根部位引入温度较低的冷气膜，从而将壁面与高温燃气隔离，对叶片表面起到冷却降温的作用。近年来，随着冷却技术的发展，异型气膜孔由于较好的冷却效率受到广泛关注。异型孔一般是由入口段的圆柱型孔与中间及出口的扩散型孔构成，研究表明异型孔的几何结构对涡轮叶片表面冷却效率有显著的影响，因此对异型孔几何结构进行测量与评定，具有重要的理论意义和实际应用价值。而目前缺乏异型孔检测与评定标准，限制了异型气膜孔的研究与应用。

发明内容

为了弥补国内缺乏异型孔检测与评定标准，本发明提供了一种基于工业CT扫描的异型气膜孔几何特征检测与评定方法。利用工业CT扫描得到异型孔试件三维点云坐标，通过剔除非异型孔轮廓点、坐标变化、截面积分、三维配准等运算处理，分别提取获得异型孔几何尺寸、形状与轮廓度等评定参数，测量精度高，算法实现相对简单，弥补了异型孔检测与评定标准的缺失，可以提高异型孔几何精度检测与评定的效率与精度。

一种基于工业CT扫描的异型气膜孔几何结构检测与评定方法，其特征在于步骤如下：

步骤1：利用工业CT扫描获得异型孔试件的三维点云数据，将所有点的三维坐标依次导出，存入齐次坐标矩阵M中，即：

其中，u为点个数，(x_mi,y_mi,z_mi)表示第i个点的X轴、Y轴、Z轴坐标，i＝1,…,u，矩阵M中包括异型孔轮廓点及试件周围平面噪声点，坐标点的存储顺序为依次按照X轴、Y轴、Z轴坐标值由小到大排序。

步骤2：分别选取坐标矩阵M中x坐标值最小的点、y坐标值最小的点和z坐标值最小的点为基准点，设定容差δ，利用下式对非异型孔轮廓点进行剔除：

其中，a、b、c分别为异型孔试件在X轴、Y轴、Z轴方向的尺寸；x_{m min}表示矩阵M中点的x坐标最小值，y_{m min}表示矩阵M中点的y坐标最小值，z_{m min}表示矩阵M中点的z坐标最小值，δ取值依据试件周围平面轮廓度确定，取值范围为[0.05,0.3]。

将剩余点的三维坐标存入齐次坐标矩阵N中，N即为异型孔的实测模型坐标矩阵，记剩余点个数为v，矩阵N表示为：

其中，v表示矩阵N中包含的点的个数，(x_nj,y_nj,z_nj)表示矩阵N中第j个点的X轴、Y轴、Z轴坐标，j＝1,…,v，矩阵N中点的坐标的存储顺序为依次按照X轴、Y轴、Z轴坐标值由小到大排序。

步骤3：提取矩阵N中x坐标值在[x_{n max}-l_x,x_{n max}]区间内的点坐标，构成异型孔入口圆柱段轮廓点坐标，其中，x_{n max}为矩阵N中点的x坐标最大值，l_x为设计模型中异型孔圆柱段X轴方向长度，将这些点的三维坐标依次存入圆柱段实测模型坐标矩阵U中：

其中，w表示矩阵U中包含的点的个数，(x_uk,y_uk,z_uk)表示矩阵U中第k个点的X轴、Y轴、Z轴坐标，k＝1,…,w，矩阵U中点的坐标的存储顺序为依次按照X轴、Y轴、Z轴坐标值由小到大排序。

按照点的x坐标值由小到大将矩阵U划分为K1段，所划分各段的x坐标值区间为[x_u1,x_u1+kΔ_x]，k＝1,2,…,K1，Δ_x为区间步长，取值范围为[0.01，0.03]，

每一段内的所有点构成一个椭圆截面，利用椭圆方程对每个椭圆截面中所有点的坐标进行拟合，得到各截面椭圆的中心点坐标，并依次存入圆柱段截面椭圆中心点坐标矩阵V中：

其中，(x_vp,y_vp,z_vp)表示矩阵V中第p个点的X轴、Y轴、Z轴坐标，p＝1,…,K1，矩阵V中点的坐标的存储顺序为依次按照X轴、Y轴、Z轴坐标值由小到大排序。

再利用最小二乘法对坐标矩阵V中所有K1个点的坐标进行直线拟合，得到异型孔中轴线的参数方程L₀：

其中，m⁰、n⁰、p⁰分别为中轴线在X轴、Y轴、Z轴方向的方向向量，(x⁰，y⁰，z⁰)为中轴线上的任意一点的点坐标，t为直线参数。

步骤4：对异型孔实测模型坐标矩阵N进行旋转坐标变换一，使异型孔中轴线与X轴平行：

S₁＝N×T₁ (7)

其中，T₁为旋转坐标变换矩阵一，S₁为进行旋转坐标变换一后的实测模型坐标矩阵。

然后，按照点的x坐标值由小到大将坐标矩阵S₁中点的坐标划分为K2段，所划分各段的x坐标值区间为[x_s1,x_s1+k′Δ′_x]，x_s1为矩阵S₁中点的x坐标最小值，k′＝1,2,…,K2，Δ′_x为区间步长，取值范围为[0.01，0.03]，

x_smax为矩阵S₁中点的x坐标最大值，每一段内的所有点即构成异型孔一个分截面，利用截面上所有点的坐标对每个分截面进行面积积分，得到其对应的分截面面积A_i，i＝1,…,K2；取所有分截面面积最小值，即得到异型孔最小截面积A_min；利用两个直线段对所有分截面面积值构成的曲线进行近似，两个直线段的交点即为异型孔圆柱段与扩散段的相交点，两个直线段在X轴的投影长度即分别为圆柱段和扩散段的孔长，记为l₁和l₂，则异型孔圆柱/扩散长度比为l₁/l₂。

步骤5：对异型孔的实测模型坐标矩阵N进行旋转坐标变换二，使异型孔中轴线与Z轴重合：

S₂＝N×T₂ (8)

其中，T₂为旋转坐标变换矩阵二，S₂为进行旋转坐标变换二后的实测模型坐标矩阵。

对坐标矩阵S₂分别进行投影坐标变换一和投影坐标变换二，使异型孔三维轮廓分别向XOZ平面和YOZ平面进行投影，得到异型孔实测模型投影坐标矩阵Q₁和Q₂：

Q₁＝S₂×F₁ (9)

Q₂＝S₂×F₂ (10)

其中，F₁为投影坐标变换矩阵一，F₂为投影坐标变换矩阵二，Q₁为进行投影坐标变换一后的实测模型投影坐标矩阵，Q₂为进行投影坐标变换二后的实测模型投影坐标矩阵。

提取投影坐标矩阵Q₁中z坐标值在[z_{q max}-l_z1,z_{q max}]区间内且x坐标值最小的坐标点，其中，z_{q max}为坐标矩阵Q₁中z坐标的最大值，l_z1为设计模型中出口扩散段轴向投影长度。利用最小二乘法对这些点的坐标进行直线拟合，得到后倾线的参数方程L₁：

其中，m¹、n¹、p¹分别为后倾线在X轴、Y轴、Z轴方向的方向向量，(x¹，y¹，z¹)为后倾线上的任意一点的点坐标。

提取投影坐标矩阵Q₂中z坐标值在[z′_{q max}-l_z1-l_z2,z′_{q max}-l_z1]区间内且y坐标值最小的坐标点，其中，z′_{q max}为坐标矩阵Q₂中z坐标的最大值，l_z2为设计模型中中间扩散段轴向投影长度，利用最小二乘法对这些点的坐标进行直线拟合，得到扩散线一的参数方程L₂：

其中，m²、n²、p²分别为扩散线一在X轴、Y轴、Z轴方向的方向向量，(x²，y²，z²)为扩散线一上的任意一点的点坐标。

提取投影坐标矩阵Q₂中z坐标值在[z′_{q max}-l_z1-l_z2,z′_{q max}-l_z1]区间内且y坐标值最大的坐标点，利用最小二乘法对这些点的坐标进行直线拟合，得到扩散线二的参数方程L₃：

其中，m³、n³、p³分别为扩散线二在X轴、Y轴、Z轴方向的方向向量，(x³，y³，z³)为扩散线二上的任意一点的点坐标。

中轴线L₀与XOY平面的夹角即为入射角α，后倾线L₁与坐标变换二后异型孔中轴线的夹角即为后倾角β，扩散线一L₂与坐标变换二后异型孔中轴线的夹角即为左扩散角γ₁，扩散线二L₃与坐标变换二后异型孔中轴线的夹角即为右扩散角γ₂。

步骤6：对齐次坐标矩阵N进行旋转平移坐标变换三，使异型孔的实测模型与设计模型的入口平面椭圆圆心及中轴线均重合，实现设计模型与实测模型坐标矩阵的粗配准。即：

S₃＝N×T₃ (14)

其中，T₃为旋转平移坐标变换矩阵三，S₃为旋转平移坐标变换三后的实测模型坐标矩阵。

步骤7：根据ICP算法对旋转平移坐标变换三后的实测模型坐标矩阵S₃和设计模型坐标矩阵S进行精配准，然后，分别提取精配准后两个模型中的异型孔入口截面、出口截面和两个相交截面的轮廓点坐标，再利用最小二乘法分别计算精配准后两个模型中这四个截面的平均偏差，即分别为入口截面、出口截面和两个相交截面的轮廓度误差。

本发明的原理及有益效果是：通过对工业CT扫描得到异型孔三维点云坐标进行非轮廓点剔除、坐标变换、截面积分、三维配准等提取得到异型孔的各分截面面积A_i最小截面面积A_min等尺寸精度表征参数，入射角α、后倾角β、左右扩散角γ₁和γ₂等形状精度表征参数，以及实测模型与设计模型的入口平面、出口平面及两个相交截面的轮廓度误差Δr等表征参数，实现了对异型孔几何尺寸、形状与轮廓度等特征的全面检测与评定。

附图说明

图1是本发明的基于工业CT扫描的异型气膜孔几何结构检测与评定方法流程图；

图2是异型(扩散型)气膜孔结构示意图及主要几何结构特征表征参数；

图中，1-入口圆柱型孔，2-出口扩散型孔，A-入口平面，B-相交平面B，C-相交平面C，D-出口平面，α-入射角，β-后倾角，γ₁-左扩散角，γ₁-右扩散角，L₀-中轴线，L₁-后倾线，L₂-扩散线一，L₃-扩散线二，l₁-异型孔圆柱段长度，l₂-异型孔扩散段长度，l_z1-异型孔出口扩散段轴向投影长度，l_z2-异型孔中间扩散段轴向投影长度，Δr_D-异型孔出口截面轮廓度误差；

图3是异型孔几何结构特征检测与评定实施例；

图中，(a)为工业CT扫描得到的三维点云图，(b)为剔除非轮廓点后的异型孔三维轮廓图，(c)为提取异型孔中轴线示意图，(d)为异型孔几何尺寸特征评定示意图，(e)为异型孔几何形状特征评定示意图，(f)为异型孔轮廓度误差评定示意图；a-异型孔试件在X轴方向尺寸，b-异型孔试件在Y轴方向尺寸，c-异型孔试件在Z轴方向尺寸，

方向基准点，P_j提取截面椭圆中心点，A_i-分截面面积，A_i/4πd²-分截面面积与设计模型圆柱段截面面积比，d-设计模型圆柱段直径。

具体实施方式

下面结合附图和实施例对本发明进一步说明，本发明包括但不仅限于下述实施例。

如图1所示，本发明提供了一种基于工业CT扫描的异型孔几何特征检测与评定方法，以超快激光加工异型气膜孔几何特征的检测与评定为对象，其中，异型气膜孔结构如图2所示。本方法通过对工业CT扫描得到异型孔三维点云坐标进行剔除非异型孔轮廓点、坐标变换、截面积分、三维配准等方法提取获得异型孔分截面面积A_i，最小截面积A_min等尺寸精度表征参数，入射角α、后倾角β及左右扩散角γ₁和γ₂等形状精度表征参数以及入口平面A、出口平面D及相交截面B和C的轮廓度误差Δr等表征参数，对异型孔几何特征进行系统评定。具体步骤如下：

步骤1：工业CT扫描获得异型孔试件的三维点云，如图3(a)所示。将三维点云导出即可获取到异型孔试件三维轮廓点坐标(X,Y,Z)，存入齐次坐标矩阵M_514920×4中，其中，514920表示扫描点个数为514920个。

其中，坐标点的存储顺序为依次按照X轴、Y轴、Z轴坐标值由小到大排序，坐标点单位为毫米(mm)。

步骤2：剔除非异型孔轮廓点坐标：由于扫描得到三维点云包括异型孔轮廓点及试件周边平面上的噪声点(即非异型孔轮廓点)，利用扫描点坐标值对非异型孔轮廓点进行剔除，分别选取坐标矩阵M_514920×4中X、Y、Z三个方向对应x、y、z坐标值最小的坐标点

(1.0611，1.4322，7.7173)，

(2.1371，1.2763，8.0290)，

(1.0686，1.4099，7.5392)作为基准点，利用下列条件进行筛选：

其中，试件在X,Y,Z方向的尺寸a，b，c经测量分别为1.41mm，1.80mm，4.12mm，如图3(a)所示。剔除点设定容差δ根据试片平面轮廓度进行选取，本实施例中试片表面经过抛光处理，平面轮廓度误差较小，取剔除点的容差δ为0.05mm即可满足要求。剔除完周边平面噪声点后，剩余点即为异型孔实测模型轮廓点，如图3(b)所示。将异型孔实测模型轮廓坐标点存入齐次坐标矩阵N_48892×4中，即N_48892×4为异型孔的实测模型坐标矩阵，48892表示异型孔实测模型轮廓点个数为48892个。

步骤3：提取拟合异型孔中轴线：异型孔由入口段圆柱型孔与出口段扩散型孔组成，本发明通过提取入口段圆柱型孔轴线来代替异型孔中轴线。首先，提取矩阵N中x坐标值在[2.1300-l_x,2.1300]区间内的点坐标，构成异型孔入口圆柱段轮廓点坐标，其中，2.1300为矩阵N中x坐标最大值，l_x为设计模型中异型孔圆柱段X向长度，本实施例中为0.6mm。将提取坐标点存入圆柱段实测模型坐标矩阵U_17641×4中，本实施例中为17641表示圆柱段轮廓点坐标个数为17641个。即：

在此基础上，按照x坐标由小到大且满足区间[1.53000,1.53000+kΔ_x],k＝1,2,…58，将坐标矩阵U划分为58段，即将异型孔圆柱段分为58个截面，Δ_x＝0.01mm为沿X向划分截面的截面间距。利用椭圆方程对单个椭圆截面坐标点矩阵中的坐标进行拟合，获得58个截面椭圆的中心点坐标矩阵V：

再利用最小二乘法对58个截面椭圆中心点坐标进行直线拟合，即可得到圆柱型孔轴线即异型孔中轴线的参数方程L₀：

其中，(0.5758，0.0133，0.8175)为中轴线的方向向量，(1.8241，2.3870，8.2984)为中轴线上的点的坐标，提取过程可见图3(c)提取异型孔中轴线示意图。

步骤4：在获得异型孔中轴线后，分别对异型孔轮廓坐标进行坐标与映射变换，并分别对异型孔几何尺寸特征、几何形状特征和截面轮廓度误差进行评定。

首先，对异型孔几何尺寸特征进行评定。对异型孔实测模型坐标矩阵N进行旋转坐标变换一，使异型孔中轴线与X轴平行，如图3(d)上图所示。

S₁＝N×T₁ (17)

其中，T₁为旋转坐标变换矩阵一，S₁为旋转坐标变换一后的实测模型坐标矩阵。根据图形变换方法，T₁为：

然后，按照步骤3方法，按照X坐标由小到大以区间长度为0.01，将坐标矩阵S₁根据区间划分为168段，即将异型孔按着截面划分为168层。根据截面坐标点对每个截面进行面积积分，得到其分截面面积A_i，i＝1,…,168，取所有分截面面积最小值，即得到异型孔最小截面积A_min＝0.0776mm²；以x坐标值为横轴、对应分截面面积值为纵轴，得到一条分截面面积的变化曲线，利用两个直线段对这条曲线进行近似，两个直线段的交点即为异型孔的圆柱段与扩散段的相交点，靠近横轴原点的直线段在横轴的投影长度记为l₁，其对应圆柱段的孔长，另外一条远离横轴原点的直线段在横轴的投影长度记为l₂，对应扩散段的孔长，异型孔圆柱/扩散长度比即为l₁/l₂＝1.08。异型孔几何尺寸特征评定过程见图3(d)异型孔几何尺寸特征评定示意图。

步骤5：对异型孔几何形状特征评定。对异型孔的实测模型坐标矩阵N进行旋转坐标变换二，使异型孔中轴线与Z轴重合，并且后倾线在XOZ平面，将异型孔“摆正”。如图3(e)中图所示。

S₂＝N×T₂ (18)

其中，T₂为旋转坐标变换矩阵二，S₂为旋转坐标变换二后的实测模型坐标矩阵；根据图形变换方法，可得旋转坐标变换矩阵T₂为：

将坐标矩阵S₂进行投影坐标变换一，将异型孔三维轮廓向XOZ平面进行投影，得到异型孔实测模型投影坐标矩阵Q₁，将坐标矩阵S₂进行投影坐标变换二，将异型孔三维轮廓向YOZ平面进行投影，得到异型孔实测模型投影坐标矩阵Q₂：

其中，F₁、F₂分别为向XOZ、YOZ平面投影矩阵，F₁，F₂分别为

提取投影坐标矩阵Q₁中z坐标满足区间[z_{q max}-0.87，z_{q max}]且x坐标值最小的坐标点，即提取XOZ平面投影轮廓上最右侧后倾段的边界点，z_{q max}为坐标矩阵Q₁中z坐标的最大值，0.87为设计模型中出口扩散段轴向投影长度。利用最小二乘法对这些坐标点进行直线拟合，得到后倾线方程参数L₁：

提取投影坐标矩阵Q₂中z坐标满足区间[z′_{q max}-1.85，z′_{q max}-0.87]且y坐标值最小/最大的坐标点，其中，z′_{q max}为坐标矩阵Q₂中z坐标的最大值，1.85为设计模型中扩散段(中间+出口扩散段)的轴向投影长度，0.87为设计模型中出口扩散段轴向投影长度。即提取YOZ平面投影轮廓上中间扩散段的左右两侧的边界点，再利用最小二乘法分别对坐标点进行直线拟合，得到扩散线一、扩散线二的参数方程L₂、L₃：

其中，(0.1956,0,0.9807)为直线L₁的方向向量，(0.2679,0,0.2244)为直线L₁上的一点的点坐标。(0,-0.1702,0.9854)为直线L₂的方向向量，(0,-0.2824,1.4964)为直线L₂上的一点的点坐标。(0,0.1529,0.9882)为直线L₃的方向向量，(0,0.2465,1.4921)为直线L₃上的一点的点坐标。

中轴线L₀与XOY平面的夹角即为入射角α，后倾线L₁与进行旋转坐标变换二后的异型孔中轴线的夹角即为后倾角β，扩散线一L₂与进行旋转坐标变换二后的异型孔中轴线的夹角即为左扩散角γ₁，扩散线二L₃与进行旋转坐标变换二后的异型孔中轴线的夹角即为右扩散角γ₂，所有异型孔形状特征参数计算结果如表1所示。

表1

入射角α	后倾角β	左扩散角γ<sub>1</sub>	右扩散角γ<sub>2</sub>
				35.15°	11.31°	8.44°	6.34°

步骤6：对异型孔轮廓度误差评定。分别利用平移变换和旋转变换组合对异型孔实测模型坐标矩阵N进行坐标变换三，使异型孔实测模型和设计模型的入口平面椭圆圆心及中轴线均重合，实现设计模型与实测模型粗配准：

S₃＝N×T₃ (23)

其中，S₃为旋转平移坐标变换三后异型孔实测模型坐标矩阵。T₃为旋转平移坐标变换三矩阵，包括平移矩阵Y和旋转矩阵R，平移距离为异型孔设计模型与实测模型两入口平面中心点坐标差，旋转角度为设计模型与实测模型两中轴线角度差。其中：

变换后实测模型与设计模型大致重合，实现粗配准。对于设计模型，如果仅有点的坐标，可以先利用等参数法提取设计模型中的异型孔轮廓点，然后按照步骤3的方法计算其中轴线及入口平面椭圆圆心点坐标。

步骤7：粗配准后，再利用ICP算法对异型孔实测模型和设计模型进行精配准，进一步提高配准精度，如图3(f)上图所示。然后，分别提取实测模型和设计模型中的异型孔入口截面A、出口截面D、相交截面B与C的轮廓点坐标，分别利用最小二乘法计算四个截面的设计与实测模型的平均偏差作为轮廓度误差Δr，其计算结果如表2所示，其中，负号表示负公差。

表2

入口轮廓度Δr<sub>A</sub>	出口轮廓度Δr<sub>D</sub>	截面B轮廓度Δr<sub>B</sub>	截面C轮廓度Δr<sub>C</sub>
				-0.0673mm	-0.0327mm	-0.0214mm	-0.0435mm

Claims (1)

1.一种基于工业CT扫描的异型气膜孔几何结构检测与评定方法，其特征在于步骤如下：

步骤1：利用工业CT扫描获得异型孔试件的三维点云数据，将所有点的三维坐标依次导出，存入齐次坐标矩阵M中，即：

其中，u为点个数，(x_mi，y_mi，z_mi)表示第i个点的X轴、Y轴、Z轴坐标，i＝1，...，u，矩阵M中包括异型孔轮廓点及试件周围平面噪声点，坐标点的存储顺序为依次按照X轴、Y轴、Z轴坐标值由小到大排序；

步骤2：分别选取坐标矩阵M中x坐标值最小的点、y坐标值最小的点和z坐标值最小的点为基准点，设定容差δ，利用下式对非异型孔轮廓点进行剔除：

其中，a、b、c分别为异型孔试件在X轴、Y轴、Z轴方向的尺寸；x_{m min}表示矩阵M中点的x坐标最小值，y_{m min}表示矩阵M中点的y坐标最小值，z_{m min}表示矩阵M中点的z坐标最小值，δ取值依据试件周围平面轮廓度确定，取值范围为[0.05，0.3]；

将剩余点的三维坐标存入齐次坐标矩阵N中，N即为异型孔的实测模型坐标矩阵，记剩余点个数为v，矩阵N表示为：

其中，v表示矩阵N中包含的点的个数，(x_nj，y_nj，z_nj)表示矩阵N中第j个点的X轴、Y轴、Z轴坐标，j＝1，...，v，矩阵N中点的坐标的存储顺序为依次按照X轴、Y轴、Z轴坐标值由小到大排序；

步骤3：提取矩阵N中x坐标值在[x_{n max}-l_x，x_{n max}]区间内的点坐标，构成异型孔入口圆柱段轮廓点坐标，其中，x_{n max}为矩阵N中点的x坐标最大值，l_x为设计模型中异型孔圆柱段X轴方向长度，将这些点的三维坐标依次存入圆柱段实测模型坐标矩阵U中：

其中，w表示矩阵U中包含的点的个数，(x_uk，y_uk，z_uk)表示矩阵U中第k个点的X轴、Y轴、Z轴坐标，k＝1，...，w，矩阵U中点的坐标的存储顺序为依次按照X轴、Y轴、Z轴坐标值由小到大排序；

按照点的x坐标值由小到大将矩阵U划分为K1段，所划分各段的x坐标值区间为[x_u1，x_u1+kΔ_x]，k＝1，2，...，K1，Δ_x为区间步长，取值范围为[0.01，0.03]，

每一段内的所有点构成一个椭圆截面，利用椭圆方程对每个椭圆截面中所有点的坐标进行拟合，得到各截面椭圆的中心点坐标，并依次存入圆柱段截面椭圆中心点坐标矩阵V中：

其中，(x_vp，y_vp，z_vp)表示矩阵V中第p个点的X轴、Y轴、Z轴坐标，p＝1，...，K1，矩阵V中点的坐标的存储顺序为依次按照X轴、Y轴、Z轴坐标值由小到大排序；

再利用最小二乘法对坐标矩阵V中所有K1个点的坐标进行直线拟合，得到异型孔中轴线的参数方程L₀：

其中，m⁰、n⁰、p⁰分别为中轴线在X轴、Y轴、Z轴方向的方向向量，(x⁰，y⁰，z⁰)为中轴线上的任意一点的点坐标，t为直线参数；

步骤4：对异型孔实测模型坐标矩阵N进行旋转坐标变换一，使异型孔中轴线与X轴平行：

S₁＝N×T₁ (7)

其中，T₁为旋转坐标变换矩阵一，S₁为进行旋转坐标变换一后的实测模型坐标矩阵；

然后，按照点的x坐标值由小到大将坐标矩阵S₁中点的坐标划分为K2段，所划分各段的x坐标值区间为[x_s1，x_s1+k′Δ′_x]，x_s1为矩阵S₁中点的X坐标最小值，k′＝1，2，...，K2，Δ′_x为区间步长，取值范围为[0.01，0.03]，

x_{s max}为矩阵S₁中点的x坐标最大值，每一段内的所有点即构成异型孔一个分截面，利用截面上所有点的坐标对每个分截面进行面积积分，得到其对应的分截面面积A_i，i＝1，...，K2；取所有分截面面积最小值，即得到异型孔最小截面积A_min；利用两个直线段对所有分截面面积值构成的曲线进行近似，两个直线段的交点即为异型孔圆柱段与扩散段的相交点，两个直线段在X轴的投影长度即分别为圆柱段和扩散段的孔长，记为l₁和l₂，则异型孔圆柱/扩散长度比为l₁/l₂；

步骤5：对异型孔的实测模型坐标矩阵N进行旋转坐标变换二，使异型孔中轴线与Z轴重合：

S₂＝N×T₂ (8)

其中，T₂为旋转坐标变换矩阵二，S₂为进行旋转坐标变换二后的实测模型坐标矩阵；

对坐标矩阵S₂分别进行投影坐标变换一和投影坐标变换二，使异型孔三维轮廓分别向XOZ平面和YOZ平面进行投影，得到异型孔实测模型投影坐标矩阵Q₁和Q₂：

Q₁＝S₂×F₁ (9)

Q₂＝S₂×F₂ (10)

其中，F₁为投影坐标变换矩阵一，F₂为投影坐标变换矩阵二，Q₁为进行投影坐标变换一后的实测模型投影坐标矩阵，Q₂为进行投影坐标变换二后的实测模型投影坐标矩阵；

提取投影坐标矩阵Q₁中z坐标值在[z_{q max}-l_z1，z_{q max}]区间内且x坐标值最小的坐标点，其中，z_{q max}为坐标矩阵Q₁中z坐标的最大值，l_z1为设计模型中出口扩散段轴向投影长度；利用最小二乘法对这些点的坐标进行直线拟合，得到后倾线的参数方程L₁：

其中，m¹、n¹、p¹分别为后倾线在X轴、Y轴、Z轴方向的方向向量，(x¹，y¹，z¹)为后倾线上的任意一点的点坐标；

提取投影坐标矩阵Q₂中z坐标值在[z′_{q max}-l_z1-l_z2，z′_{q max}-l_z1]区间内且y坐标值最小的坐标点，其中，z′_{q max}为坐标矩阵Q₂中z坐标的最大值，l_z2为设计模型中中间扩散段轴向投影长度，利用最小二乘法对这些点的坐标进行直线拟合，得到扩散线一的参数方程L₂：

其中，m²、n²、p²分别为扩散线一在X轴、Y轴、Z轴方向的方向向量，(x²，y²，z²)为扩散线一上的任意一点的点坐标；

提取投影坐标矩阵Q₂中z坐标值在[z′_{q max}-l_z1-l_z2，z′_{q max}-l_z1]区间内且y坐标值最大的坐标点，利用最小二乘法对这些点的坐标进行直线拟合，得到扩散线二的参数方程L₃：

其中，m³、n³、p³分别为扩散线二在X轴、Y轴、Z轴方向的方向向量，(x³，y³，z³)为扩散线二上的任意一点的点坐标；

中轴线L₀与XOY平面的夹角即为入射角α，后倾线L₁与坐标变换二后异型孔中轴线的夹角即为后倾角β，扩散线一L₂与坐标变换二后异型孔中轴线的夹角即为左扩散角γ₁，扩散线二L₃与坐标变换二后异型孔中轴线的夹角即为右扩散角γ₂；

步骤6：对齐次坐标矩阵N进行旋转平移坐标变换三，使异型孔的实测模型与设计模型的入口平面椭圆圆心及中轴线均重合，实现设计模型与实测模型坐标矩阵的粗配准，即：

S₃＝N×T₃ (14)

其中，T₃为旋转平移坐标变换矩阵三，S₃为旋转平移坐标变换三后的实测模型坐标矩阵；

步骤7：根据ICP算法对旋转平移坐标变换三后的实测模型坐标矩阵S₃和设计模型坐标矩阵S进行精配准，然后，分别提取精配准后两个模型中的异型孔入口截面、出口截面和两个相交截面的轮廓点坐标，再利用最小二乘法分别计算精配准后两个模型中这四个截面的平均偏差，即分别为入口截面、出口截面和两个相交截面的轮廓度误差。

Images