基于射线衰减能量场的无损检测缺陷提取、识别方法

摘要

基于射线衰减能量场的无损检测缺陷提取、识别方法，首先，根据射线检测以及数字成像建立底片图像灰度与衰减能量场分布之间的关系，从射线底片图像中获取缺陷边界；然后，基于射线衰减能量场对所述缺陷提取方法得到的缺陷进行类型识别；最后，根据能量衰减原理，对所述缺陷类型识别方法得到的缺陷进行体积测量。由于本发明利用射线能量穿透不同材质所表现的衰减特征，建立射线能量衰减特征与底片图像灰度的对应规律，利用计算机对灰度图像的解析度远远超过人眼的优势，从而实现产品零部件缺陷的自动提取、识别以及三维体积测量。

说明书

技术领域

本发明涉及一种数字图像处理中基于射线衰减能量场的无损检测缺陷提取、识别方法。

背景技术

利用工业射线底片图像判定缺陷主要包括缺陷边界提取和缺陷类型识别两个方面。由于零部件内部缺陷的多样性与复杂性，目前，大部分缺陷类型识别主要是以缺陷所形成的影像底片作为判别依据，依靠专业人员的实际经验进行判定，造成无损检测技术人员工作量大，工作强度高，错判漏判现象时有发生。此外，无损射线检测行业涉及到探伤工艺、铸造工艺、材料成形工艺、焊接工艺以及表面处理工艺，培养一名合格的评片人员大约需要三年时间。因此，利用专业装置和计算机进行辅助缺陷识别，减轻工作人员劳动强度，提高缺陷识别率成为必须。目前利用射线底片图像进行零部件缺陷识别主要存在以下问题：

缺陷提取方面，目前主要有边界方法、区域方法和聚类方法三种，边界方法主要利用图像相邻像素间的灰度突变进行微分运算，从而检测出不同区域的边界，这种方法对噪声非常敏感且无法识别微小缺陷。区域方法主要是根据像素之间的相似度，把像素组合成不同的区域，这类方法容易产生过分割，需要进一步消除。聚类方法主要是在特征空间上，把相似的点聚类在一起，然后在图像中用不同的类号来标注出不同的区域，这类方法的主要缺点是聚类所需的类总数一般都不可知，而且往往会忽视相邻像素之间的空间和视觉上的相互关系。

缺陷类型识别方面，目前主要采用神经网络和专家系统，利用神经网络识别缺陷首先需要提取表征缺陷类型的特征，如缺陷形状、大小、位置等，然后利用样本对神经网络进行训练，用训练好的神经网络对未知缺陷进行识别。这种方法需要大量的缺陷样本，训练时间也比较长，识别的准确度与样本的数量以及特征的选取有很大关系，目前，国外资料显示神经网络对气孔的识别率为92.39％，对于其他缺陷的识别率相对较低。在利用专家系统识别方面，由于知识是决定一个专家系统性能是否优越的关键因素，专家系统要达到领域专家的水平，就必须掌握领域专家处理问题时所使用的大量专门知识，特别是经验知识，从专家对大量实例的分析中获取专家解决问题的思路、知识、经验及规则。因此，获取知识并把知识表达成专家系统可用的形式，是专家系统开发中的主要瓶颈之一。建立一个实用的专家系统需要一个相当长的时间。更为重要的是，目前采取的所方法，无论是人工识别还是计算机判别(神经网络、专家系统)，其共同点都是通过分析缺陷的边界形状特征来识别缺陷类型，很少考虑缺陷内部区域灰度分布，包含在缺陷区域内部与缺陷类型密切相关的大量宝贵信息被浪费。

在缺陷测量方面，二维测量如缺陷长度和面积测量已经十分成熟，但利用图像灰度实现缺陷三维体积测量基本是一个空白。

发明内容

本发明的目的是为了克服上述现有技术仅仅利用与缺陷类型部分相关的缺陷边界特征而忽略了包含在缺陷内部的与缺陷类型更为密切的灰度信息的缺点，提出了一种基于射线衰减能量场的无损检测缺陷提取、识别方法，该方法利用衰减能量场与图像灰度分布之间的对应关系，利用计算机充分挖掘包含在缺陷内部的与缺陷类型更为密切的灰度特征信息，实现缺陷自动提取、识别以及三维体积测量。

为达到上述目的，本发明采用技术方案是：

1)依据射线穿过不同物质的能量衰减规律，建立射线检测底片数字化图像灰度与射线衰减能量场的关系，确定射线底片图像中的缺陷边界；

2)基于射线能量场的衰减特征对所确定的缺陷进行类型识别；

3)根据射线能量衰减的强度以及与图像灰度的对应关系，对缺陷体积进行测量。

所述第一步骤包括：

1)缺陷种子列搜索步骤：利用种子列判别条件识别定位缺陷，即通过底片图像列灰度曲线即由某列像素灰度值构成的曲线波谷的宽度、深度以及与波谷相邻两波峰的高度差，确定缺陷波谷，从而识别定位缺陷；

2)缺陷概略分割步骤：据缺陷边缘的连续性，限定每列的搜索范围，分别向种子列两侧搜索缺陷；

3)缺陷细节提取步骤：通过概略分割结果重组，将分布于相邻几个不同区域的同一缺陷的几个不同部分合并成一个区域，利用直方图统计，对所述区域逐层进行直方图统计，将各层出现频率最高的灰度值绘成灰度曲线，根据所述灰度曲线确定缺陷分割阈值，依据阈值对所述区域进行二次分割，通过边缘完整性检验，确定缺陷边界是否存在缺口，如果存在，依据所述阈值，沿缺口连线法线方向逐层搜索边界，从而获取完整的缺陷边界；

所述第二步骤包括：

1)缺陷材质类型判别步骤：通过实验建立三类不同材质：金属、氧化物、气体在不同射线强度下的灰度——材质库，找出相同射线强度下缺陷平均灰度值对应的T_A、T_B、T_C，通过缺陷比较当量尺寸和T_A、T_B、T_C，确定缺陷材质类型；

2)缺陷特征提取步骤：对所述缺陷细节提取获取缺陷区域灰度曲面进行三维形貌特征和二维形状位置特征提取，三维形貌特征包括：单位长度交叉点数、深度宽度比、曲面凸凹程度、内部区域小波谷数量、小波谷位置、小波谷长度，二维形状位置特征包括：缺陷长宽比、缺陷重心坐标；

3)缺陷形貌特征库提炼步骤：通过对样本缺陷灰度曲面的三维形貌特征和二维形状位置特征进行提取，分析、建立缺陷形貌特征库；

4)缺陷类型识别步骤：对缺陷区域灰度曲面的三维形貌特征和二维形状位置特征进行提取，结合缺陷材质类型，依据缺陷形貌特征库对提取特征进行分析，实现缺陷类型自动识别；

所述第三步骤包括：

缺陷体积测量步骤：利用缺陷内部各点与缺陷边界的灰度差以及由缺陷材质和透照条件决定的射线衰减系数，根据公式ΔT＝k·Δv/(μ-μ′)得出缺陷内部各点在射线透照方向上的尺寸，从而得到缺陷的体积。

利用射线能量穿透不同材质所表现的衰减特征，建立射线能量衰减特征与底片图像灰度的对应规律，利用计算机对灰度图像的解析度远远超过人眼的优势，对产品零部件射线无损检测缺陷进行自动提取、识别以及三维体积测量。

附图说明

图1是射线检测基本原理；

图2(a)是实验图像，其中：x为图像的行坐标，y为图像的列坐标；

图2(b)是图2(a)的灰度曲面图，其中x为图2(a)的行坐标，y为图2(a)的列坐标，v代表像素点灰度值；

图3是缺陷种子列示意图，其中：i表示种子的位置，x_t^ij，x_p1^ij，x_p2^ij分别表示种子列(图像的第i列)上第j个缺陷的波谷及左、右两相邻波峰的位置坐标，v_t^ij，v_p1^ij，v_p2^ij表示图像上x_t^ij，x_p1^ij，x_p2^ij对应点像素的灰度值；

图4是相邻列缺陷边界点位置示意图，其中：x为行坐标，y为列坐标；

图5是概略分割流程图；

图6是缺陷概略分割实验结果；

图7(a)是对缺陷区域逐层进行灰度直方图统计示意图；

图7(b)是由图7(a)各层出现频率最高的灰度值构成的灰度曲线图；

图8是缺陷细节提取流程图；

图9(a)是本发明对图2(a)的缺陷提取结果；

图9(b)是Canny方法对图2(a)的缺陷提取结果；

图9(c)是分水岭方法对图2(a)的缺陷提取结果；

图9(d)是水平集方法对图2(a)的缺陷提取结果；

图10(a)是单位长度交叉点示意图。

图10(b)是深宽比示意图。

图10(c)曲面凸凹程度示意图。

图10(d)是小波谷示意图。

图11是缺陷形貌特征库创建流程图。

图12是缺陷灰度——材质曲线示意图。

图13是缺陷类型识别流程图。

图14是缺陷像素点示意图。

具体实施方式

下面结合附图对本发明作进一步详细说明。

图1是射线检测基本原理图，图中T、ΔT分别表示工件基体和缺陷在射线透过方向上的尺寸，I₀表示入射射线强度，I₁、I₂分别表示穿透工件基体和缺陷部位的射线强度，μ、μ′分别表示工件基体和缺陷对射线的衰减系数。射线穿过工件后经过衰减的能量场(射线衰减能量场)分布如下：

在正常曝光情况下底片黑度值随射线照射量对数的增加而呈线性增长，底片经CCD数字化处理后，得到的数字化图像灰度值与原底片灰度值呈线性关系。由此得出底片图像灰度与衰减能量场关系如下：

v＝k·lg(I)+m    (2)

其中：k，m为常数且k＜0，可以看出图像的灰度值随射线衰减能量强度对数的增加而呈线性减小。

图2(a)是实验图像，以焊缝γ射线底片局部数字图像(85像素×58像素)为例对本发明方法进行验证。

图2(b)是图2(a)对应的灰度曲面，灰度曲面由图像像素点坐标x，y及对应灰度值v三个坐标构成，由式(2)可知图像各像素点的灰度值与对应的衰减能量场射线强度具有精确的映射关系，基于此，灰度曲面的表面状况可以直观的反映射线衰减能量场的分布。在同等射线强度下除了“夹钨”缺陷(钨的射线衰减系数比基体射线衰减系数要高的多，夹钨缺陷对应的灰度曲面为凸台)外缺陷材质的射线衰减系数比基体的射线衰减系数要小，对应灰度曲面上为凹坑，且凹坑的表面形貌特征反映了缺陷的类型，根据这些特征对进行缺陷提取、类型识别以及体积测量。

图3是缺陷种子列示意图，x_t^ij，x_p1^ij，x_p2^ij分别表示种子列(图像的第i列)上第j个缺陷的波谷及左、右两相邻波峰的位置，v_t^ij，v_p1^ij，v_p2^ij表示图像上x_t^ij，x_p1^ij，x_p2^ij对应像素点的灰度值v，缺陷种子列(图像中可以用来判断缺陷的列)需同时满足以下4个判别条件：

1)缺陷波谷具有一定的深度：

Min{Δv₁，Δv₂}＞h

其中：${\mathrm{\Delta v}}_1=v_{p2}^{\mathrm{ij}}-v_t^{\mathrm{ij}},$${\mathrm{\Delta v}}_2=v_{p1}^{\mathrm{ij}}-v_t^{\mathrm{ij}};$h：波谷与相邻两波峰灰度差阈值。

2)缺陷波谷具有一定的宽度：

Min{Δx₁，Δx₂}＞w

其中：${\mathrm{\Delta x}}_1=x_{p2}^{\mathrm{ij}}-x_t^{\mathrm{ij}},$${\mathrm{\Delta x}}_2=x^t-x_{p1}^{\mathrm{ij}};$w：波谷与相邻两波峰的距离阈值。

3)波谷到左右两相邻波峰的距离差不能超过波谷宽度的1/α：

|Δx₁-Δx₂|＜Δw

其中：$\mathrm{\Delta w}=(x_{p2}^{\mathrm{ij}}-x_{p1}^{\mathrm{ij}})/\alpha ;$α∈N，α＞1。

4)波谷左右两相邻波峰的灰度差不能超过两个波峰平均灰度的1/β：

|Δv₁-Δv₂|＜Δh

其中：$\mathrm{\Delta h}=(v_{p1}^{\mathrm{ij}}+v_{p2}^{\mathrm{ij}}-{2v}_t^{\mathrm{ij}})/\beta :\beta \in N,\beta >1.$

图4是相邻列缺陷边界点位置示意图，图中第i列第j个缺陷的上下边界点坐标分别为x_p1^ij，x_p2^ij，根据缺陷边界的连线性，i-1、i+1列上缺陷的上下边界点坐标应在区间[x_p1^ij-1，x_p2^ij+1]内，考虑到噪声等因素的影响，缺陷在i-1、i+1列上的边界点坐标x^i-1j，x^i+1j应满足如下条件：

$x^{i-1j},x^{i+1j}\in [x_{p1}^{\mathrm{ij}}-\delta ,x_{p2}^{\mathrm{ij}}+\delta ]---\left(3\right)$

其中；δ是偏移量，取值在3～5之间。

图5是概略分割流程图，具体描述如下：

步骤S1，读入一幅图像。

步骤S2，对图像进行低通滤波。

步骤S3，依据缺陷种子列判别条件，从图像第一列开始向后逐列搜索缺陷种子列。

步骤S4，判断是否找到缺陷种子列。

上述步骤S4的判定为“是”，即找到一个缺陷种子列，进行步骤S5，记录种子列位置i，缺陷个数nDef，每个缺陷的波谷位置x_t^ij、相邻左波峰的位置x_p1^ij、相邻右波峰的位置x_p2^ij，以及x_t^ij，x_p1^ij，x_p2^ij对应像素点的灰度值v_t^ij，v_p1^ij，v_p2^ij，其中：j＝1，…，nDef。

上述步骤S4的判定为“否”，即没有找到缺陷种子列，进行步骤S17，搜索结束。

步骤S6，定位到第j个缺陷，即准备对第j个缺陷进行概略分割。

步骤S7，将i的值赋值给k，即k的值为缺陷列的位置。

步骤S8，向前搜索k-1列上满足条件(3)区间内的缺陷内部最小灰度值点和上下边界点即灰度曲线的波谷以及相邻两波峰的坐标x_t^k-1j，x_p1^k-1j，x_p2^k-1j。

步骤S9，判断以下两个条件：1)搜索没有到达图像左边界(k-1＞0)，2)k-1列的两波峰值大于k列的波谷值$(v_t^{\mathrm{kj}}<v_{p1}^{k-1j}\mathrm{and}{v}_t^{\mathrm{kj}}<v_{p2}^{k-1j}),$是否同时满足。

上述步骤S9判定为“否”，进行布骤S10，继续向前搜索下一列(k＝k-1)。

上述步骤S9判断为“是”，进行步骤S11，即种子列左侧缺陷部分搜索完毕，从种子列开始(k＝i)，搜索种子列右侧缺陷部分。

步骤S12，向后搜索k+1列上满足条件(3)区间内的缺陷内部最小灰度值点和上下边界点即灰度曲线的波谷以及相邻两波峰的坐标x_t^k+1j，x_p1^k+1j，x_p2^k+1j。

步骤S13，判断以下两个条件：1)搜索没有到达图像右边界(k＜N)，2)k+1列的两波峰值大于k列的波谷值$(v_t^{\mathrm{kj}}<v_{p1}^{k+1j}\mathrm{and}{v}_t^{\mathrm{kj}}<v_{p2}^{k+1j}),$是否同时满足。

上述步骤S13判定为“否”，进行布骤S14，继续向后搜索下一列(k＝k+1)。

上述步骤S13判断为“是”，进行步骤S15，即第j个缺陷搜索完毕，准备搜索j+1个缺陷(j＝j+1)。

步骤S16，判断种子列上的nDef个缺陷是否都已搜索完毕(j＜nDef+1)。

上述步骤S16判定为“是”，即第j各缺陷未进行搜索缺陷，进行步骤S6，开始搜索第j个缺陷。

上述步骤S16判定为“否”，即种子列i上的nDef个缺陷都已搜索完毕，进行步骤S3，从第i列开始向后搜索种子列，搜索区域不包括已进行缺陷搜索的区域。

图6是图2(a)概略分割结果，由上下两条边界构成的一个区域。

图7(a)是对缺陷区域逐层进行灰度直方图统计示意图，

k₁，k₂，…，k_j，…，k_ε，(j＝1，…，ε)为缺陷的层数。

图7(b)是由图7(a)k₁，k₂，…，k_j，…，k_ε各层出现频率最高的灰度值构成的灰度曲线，根据该灰度曲线确定缺陷的分割阈值。

图8是缺陷细节提取流程图，参照流程图对缺陷细节提取过程说明如下：

步骤S1，读入一幅图像的概略分割结果。

步骤S2，对概略分割结果进行重组，即将分布于相邻几个区域内的同一缺陷的不同部分合成一个区域。

步骤S3读取一个缺陷区域。

步骤S4，对该区域逐层进行直方图统计(图7a)，将各层出现频率最高的灰度值绘制成灰度曲线(图7b)。

步骤S5，根据灰度曲线找出最佳分割阈值(图7b)。

步骤S6，依据上述步骤S5得到的阈值对步骤S3读取的区域进行二次分割，然后对分割获取边界进行完整性检验，如果存在缺口，按照所述阈值沿缺口连线法线方向逐层搜索缺陷边界，从而获得完整的缺陷边界。

步骤S7，判断是否图像上的所有区域都已进行了细节提取。

上述步骤S7判断为“否”，进行步骤S3，即对下一个缺陷进行细节提取。

上述步骤S7判断为“是”，细节提取结束。

图9(a)是本发明对图2(a)的缺陷提取结果。

图9(b)是Canny方法对图2(a)的缺陷提取结果。

图9(c)是分水岭方法对图2(a)的缺陷提取结果。

图9(d)是水平集方法对图2(a)的缺陷提取结果。

由图9(a)、(b)、(c)、(d)的对比结果可以看出，本发明的提取效果明显好于Canny方法和分水岭方法，与水平集方法相比本发明方法的提取结果保留了更多的细节信息，且具有明显的时间优势。

图10是缺陷灰度曲面三维形貌特征示意图，三维形貌特征包括：单位长度交叉点数、深度宽度比、曲面凸凹程度、内部区域小波谷数量、小波谷位置、小波谷长度。

图10(a)是单位长度交叉点示意图，图中左侧为缺陷某一列的灰度曲线(列灰度曲线)，右侧为构成缺陷灰度曲面的所有列灰度曲线的波谷谷底连线，将该曲线与首尾两点连线的交叉点数N除以首尾两点间的距离L，得到长度交叉点个数：n＝N/L。

图10(b)是深宽比示意图，图中左侧为缺陷某一列的灰度曲线(列灰度曲线)，右侧为构成缺陷灰度曲面的所有列灰度曲线的波谷深度H与波谷宽度W的比值绘成一条曲线，直线$\frac{H}{W}=1$将曲线分成三个区域，宽度分别为：D₁、D₂、D₃，D₁、D₂、D₃三者之间的比例反映了灰度曲面的深宽比。

图10(c)曲面凸凹程度示意图。图中的d1、d2分别为左右两段灰度曲线距连线的最远距离，将所有列灰度曲线的d/D(d＝(d₁+d₂)/2，D＝(D₁+D₂)/2)求平均值，该平均值反映了灰度曲面的凸凹程度。

图10(d)是小波谷示意图。图中显示了缺陷内部区域灰度曲面的所有小波谷数量、位置，长度。

图11是缺陷形貌特征库创建流程图，参照流程图对缺陷形貌特征库创建过程说明如下：

步骤S1，读取缺陷样本。

步骤S2，提取所述步骤S1读入缺陷的三维形貌特征。

步骤S3，提取所述步骤S1读入缺陷的形状位置特征。

步骤S4，分析归纳三维形貌特征和形状位置特征。

步骤S5，根据所述步骤S4的分析归纳结果，建立缺陷形貌特征库。

图12是缺陷——灰度材质曲线示意图，图中V代表灰度，T代表射线透过方向材质的尺寸。通过实验建立三种不同材质：金属、氧化物、气体在不同射线强度下的灰度——材质库，将缺陷材质的当量尺寸与相同射线强度下的平均灰度对应的T_A、T_B、T_C进行比较，确定缺陷的材质类型。

图13是缺陷类型识别流程图，参照流程图对缺陷类型识别过程说明如下：

步骤S1，读入一幅图像的细节提取结果。

步骤S2，读取一个缺陷。

步骤S3，对上述步骤S2读取缺陷的灰度曲面进行三维形貌特征提取。

步骤S4，对上述步骤S2读取缺陷进行形状位置特征提取。

步骤S5，计算缺陷区域的平均灰度和当量尺寸(表征射线透过方向上缺陷尺寸的当量值)。

步骤S6，依据灰度——材质库，得出缺陷材质类型。

步骤S7，模拟专家识别缺陷的机理，利用缺陷形貌特征库并结合缺陷材质类型对缺陷特征进行分析，实现缺陷自动识别。

步骤S8，判断是否图像中的所有缺陷都已进行了识别。

上述步骤S8判定为“否”，进行步骤S2，即识别下一个缺陷。

上述步骤S8判定为“是”，识别过程结束。

图14是缺陷像素点示意图，图中每个方格代表一个像素点。通过底片灰度与射线衰减能量场之间的关系(式2)以及射线透过方向上缺陷尺寸与衰减能量场射线强度之间的关系(式1)得出射线透过方向上缺陷尺寸与缺陷灰度差(缺陷内部对各点与缺陷边界的灰度差)之间的关系如下：

ΔT＝k·Δv/(μ-μ′)    (4)

式中，k是常数，μ是基体材质的射线衰减系数，μ′是缺陷材质的射线衰减系数，根据式(4)计算得出缺陷内部各点射线透过方向上尺寸ΔT进而得到缺陷体积。

本发明设计了一种专门针对射线检测图像的缺陷提取、缺陷类型识别以及缺陷体积测量方法。缺陷提取方法比目前用于缺陷提取的分水岭方法、水平集方法，canny方法效果要好，图5(a)(b)(c)(d)以焊缝缺陷射线图像为例对上述算法的提取效果进行了比较。表1给出了该焊缝缺陷图像采用不同方法进行缺陷提取的平均处理时间(各个算法对图2a的25次运算耗时平均值)，可以看出本发明方法的平均处理时间与Canny方法、分水岭方差不多，但水平集方法相比，具有明显的优势。为了评估本发明的提取精度，采用平均相似性准则进行评估：

$S=\frac{1}{N}\underset{i=1}{\overset{N}{\Sigma }}\frac{|U_g\cap U_s|}{|U_g\cup U_s|}---\left(5\right)$

式中：U_g和U_s分别表示理想提取和算法提取得到的像素点集合；N是图像的个数。4个有经验的研究人员分别对10幅缺陷图像进行手动分割，将每幅图像的4个不同的分割结果进行平均，将平均结果作为每幅图像的理想缺陷提取，根据式(5)计算平均相似性，表2给出了比较结果。总体来看，本发明方法在平均处理时间和提取精度上均优于三种与之比较的方法。缺陷识别方法突破了传统缺陷识别过程中主要依靠影像或图像边界信息进行缺陷识别的局限性，利用计算机对图像的高解析度优势充分挖掘隐藏在缺陷内部灰度区域中与缺陷类型具有密切关系的灰度特征，依据缺陷形貌特征并结合缺陷材质类型对缺陷进行分行，实现缺陷类型自动识别。缺陷体积测量方法充分利用射线能量衰减原理，通过灰度差与能量差以及能量差与射线透过方向上缺陷尺寸之间的关系，得出缺陷体积。

实验的硬件环境为Celeron处理器，主频为2.80GHz，软件环境为windows XP操作系统，Visual C++.NET2003编程语言。实验图像(图2a)为国内某大型汽轮机厂高压阀焊缝γ射线底片局部数字图像(58像素×85像素)。

表1平均处理时间

算法    t/msCanny方法分水岭方法水平集方法本文方法    6    10    1670    5

表2提取精度

算法    SCanny方法分水岭方法水平集方法本文方法    0.7462    0.9160    0.9965    0.9986