技术领域
[0001] 本
发明涉及一种射线追踪式超声Lamb波缺陷层析成像方法，属于
无损检测技术领域。
背景技术
[0002] 利用超声Lamb波缺陷层析成像技术可以快速、有效地获得缺陷的轮廓和尺寸等具体信息。该技术继承了传统超声Lamb波检测的诸多优点。Lamb波走时(Time-of-flight，TOF，指Lamb波在收发换能器对之间的传播时间)跨孔层析成像技术，利用线性换能器阵列进行扇束投影和成像，是一种基于
迭代法的高效、快速的层析成像方法。传统的Lamb波TOF跨孔层析成像均是基于直射线理论，即假设Lamb波射线在传播过程中遇到缺陷后仍然沿直线传播。但是，当平板中存在强散射缺陷时，Lamb波射线传播方向会发生偏转，即Lamb波射线实际将沿弯曲路径传播。此时，已有的直射理论将不再适用，并将导致层析成像的
精度下降。上述问题是制约Lamb波走时跨孔层析成像技术发展的一个
瓶颈问题。中国
专利文献公开了一种“刚制储罐
底板的导波在线检测方法”该技术涉及一种
钢质储罐底板的导波在线检测方法，采用超声
探头阵列收发导波，通过
滤波反投影算法重建层析图像，可快速实现对储罐底板的全面检测和缺陷分析，但该技术局限于矩阵变换理论，只能判断出缺陷的存在
位置及严重程度，而对其具体尺寸和轮廓并不能达到十分精确的估计。
发明内容
[0003] 本发明的目的是提供一种射线追踪式超声Lamb波缺陷层析成像方法，该方法按如下步骤进行：
[0004] 步骤1：在待测材料上
选定成像区域，将所述成像区域划分成N1×N2个小网格；在成像区域的一侧设置M个发射EMAT(电磁
声换能器)，对侧相同位置处设置M个接收EMAT；其中，N1、N2、M为正整数；
[0005] 步骤2：激励第1～M个发射EMAT全向激发A0模态Lamb波，M个接收EMAT依次接收Lamb波；
[0006] 步骤3：利用伪Wigner-Ville分布(PWVD)对M×M个EMAT检测
波形进行时频分析和模态识别；
[0007] 步骤4：提取M×M个检测波形的时频分析结果，记录M×M个A0模态Lamb波的走时T1～TM×M；
[0008]
[0009] 步骤5：使用联合
迭代重建算法确定每个小网格的慢度及缺陷分布；
[0010] 步骤6：在成像区域中建立三维直
角坐标系，然后建立初始发射角度，利用外推法逐步对经过缺陷区域的全部射线路径进行追踪，对于射线i前进方向上的一点P(x，y)，得到Lamb波的相速度值cp和 和 其中，x、y为P点沿x方向、y方向的坐标值[0011] 步骤7：将得到的cp、 和 代入射线追踪外推公式，找到新的外推点；
[0012] 步骤8：判断新找到的外推点坐标值是否到达或超出子网格区域边界的坐标值，若否，则返回步骤6；若是，则结束对射线i的射线追踪，记录此时外推点 的坐标作为终点坐标；
[0013] 步骤9：判断终点坐标是否到达或接近接收EMAT，若是，则得到射线i的预期最短走时路径；若否，则改变初始发射角的值，返回步骤6，直至找到连接收发EMAT的射线i的最短走时路径；
[0014] 步骤10：重复步骤6～9，直至找到经过缺陷区域的全部射线的最短走时路径。
[0015] 步骤11：将得到的全部最短射线路径作为已知量代入步骤5联合迭代重建算法，再次得到网格慢度，以及慢度曲线在缺陷区域的突变情况，判断缺陷大小及分布。
[0016] 步骤12：判断缺陷分布的成像精度是否满足要求，若是，则得到预期结果，若否，则重复步骤2～11。
[0017] 收发EMAT直径为20～80mm，相邻EMAT中心间距为20～100mm。
[0018] 射频功率
放大器激励发射EMAT，射频
功率放大器的激发
频率为50～500kHz。
[0019] 激发的Lamb波为单一的A0模态Lamb波。
[0020] 联合迭代重建算法确定每个网格的慢度及缺陷分布的方法为：
[0021]
[0022] 其中，Sj为待求的第j个网格的慢度；Lij为第i条投影射线在第j个网格中的长度；Ti为第i条投影射线的实测走时；n为正整数，且n＝N1×N2；m为正整数，且m＝M×M。
[0023] 初始发射角度的范围是0～180度。
[0024] 利用外推法逐步对射线路径进行追踪的方法为：以外推点为中心获取16个
像素；得到相应的16个cp值；将cp值代入二维三次样条插值公式，对Lamb波 cp的分布情况进行二维曲面拟合，求解每一步外推所需的 和 求解公式为：
[0025]
[0026] 其中：cp是相速度，x，y分别表示16个像素中心点沿x轴方向、y轴方向的坐标值；Ckl表示与16个像素相对应的待求解系数；k,l＝0,1,2,3。
[0027] 根据射线追踪外推公式，寻找新的外推点的方法为：
[0028]
[0029] 其中，yk表示经过k步外推获得的射线i上的离散点沿y轴方向的坐标值；yk+1、yk-1分别表示上述离散点前、后两点沿y轴方向的坐标值；Δx表示射线外推的固定步长，且[0030] 改变初始发射角的值的方法为：每次增加或减少2～6度。
[0031] 本发明的有益效果是：本发明将RT迭代算法模型与TOF跨孔层析成像技术相结合，综合两者优点，通过正反问题的交替反复求解，对正问题中的射线路径加以修正，以还原缺陷的真实分布。本发明计算准确、高效、快速，解决了基于直射线理论的传统层析成像方法成像精度低的问题，具有广阔的应用前景；同时避免了基于直射线理论的传统TOF跨孔层析成像方法成像精度低的问题，计算准确、高效、快速，为提高缺陷成像精度提供重要依据。
附图说明
[0032] 图1为射线追踪式超声Lamb波层析成像
流程图。
[0033] 图2为成像区域网格划分示意图。
具体实施方式
[0034] 本发明提出的方法，其实质是建立一种全新的射线追踪(Ray Tracing，RT)迭代算法模型，并将其与TOF跨孔层析成像技术相结合，综合两者优点，通过正反问题的交替反复求解，对正问题中的射线路径加以修正，以获得更为真实的缺陷分布情况。下面结合
实施例对本发明做进一步说明：
[0035] 步骤1：取厚度为1～5mm的待测
铝板，在上表面选择一
块面积为A×B的矩形区域作为成像区域进行实验，其中，A、B为正实数，且均大于30mm。在所述区域内加工一个直径大于等于30mm、深为1～5mm的标准人造
腐蚀坑缺陷。然后将所述成像区域划分成N1×N2个小网格(即N1×N2个像素)。在成像区域的一侧等间距设置M个发射EMAT(电磁声换能器)，对侧相同位置处设置相同数量的接收EMAT。所述的EMAT直径为20～80mm，同侧相邻EMAT中心间距为20～100mm。
[0036] 步骤2：选择一台射频功率放大器用于激励EMAT，令发射端中所有M个发射EMAT均被全向激发A0模态Lamb波，所有M个接收EMAT进行接收，激发频率为50～500kHz。所述的Lamb波为单一的A0模态Lamb波，不含有A0以外其他模态Lamb波。
[0037] 步骤3：利用伪Wigner-Ville分布(PWVD)对M×M个EMAT检测波形进行时频分析和模态识别；
[0038] 步骤4：提取所有M×M个检测波形的时频分析结果，并记录M×M个A0模态Lamb波的走时，记为T1～TM×M。
[0039] 步骤5：使用联合迭代重建算法(SIRT)公式确定每个小网格的慢度(速度的倒数)及缺陷的分布：
[0040]
[0041] 其中，Sj为待求的第j个小网格的慢度；Lij为第i条投影射线在第j个小网格中的长度；Ti为第i条投影射线的实测走时，n＝N1×N2，m＝M×M；
[0042] Lij通过几何关系计算得出：已知成像区域面积，小网格边长，收发EMAT数量、位置和间距；按照直射线理论计算得出射线路径，可知射线从发射端到接收端经过的每一个小网格中的长度。
[0043] 步骤6：首先在成像区域建立一个三维直角坐标系，其中，x、y代表射线上的一点在二维平面上的位置，z代表该点的相速度cp，然后进行子网格区域划分，每个子网格区域包含4×4＝16个像素，接着建立一个初始发射角度，范围在0～180度之间；
[0044] 利用外推法逐步对经过缺陷区域的全部射线路径进行追踪；
[0045] 对于第i条射线，利用外推法逐步对其射线路径进行追踪。根据Lamb波的传播特性可知，当遇到腐蚀缺陷时，其相速度cp会随板厚的变化而变化。因此，在外推点区域，以外推点为中心，获取4×4区域的16个像素，得到对应16个A0模态Lamb波的相速度cp值，代入二维三次样条插值(2D-CSI)公式：
[0046]
[0047] 构成16个线性方程，对Lamb波cp的分布情况进行二维曲面拟合，进而精确求解每一步外推所需的 和
[0048] 其中，cp是相速度，x，y表示沿追踪方向x、y的坐标值；Ckl(k,l＝0,1,2,3)表示16个待求解系数，16个像素中每一个像素都对应一个x、y值。
[0049] 步骤7：将步骤6得到的cp、 和 代入射线追踪(Ray Tracing， RT)外推公式，完成一步外推：
[0050]
[0051] 其中，yk表示经过k步外推获得的射线i上的离散点沿y轴方向的坐标值；
[0052] yk+1、yk-1分别表示上述离散点前、后两点沿y轴方向的坐标值；
[0053] Δx表示射线外推的固定步长，具体有：
[0054]
[0055] 步骤8：判断新找到的外推点坐标值是否到达或超出子网格区域边界坐标值，若不是，则返回步骤5；若是，则结束对射线i的射线追踪，记录下终点坐标(xf，yf)。
[0056] 步骤9：判断终点坐标(xf，yf)是否到达或接近目标接受位置，若是，则得到射线i的预期最短走时路径；若不是，则将初始发射角每次增加或减少2～6度，再重复步骤6～8，直到找到连接收发EMAT的射线i的最短走时路径。
[0057] 步骤10：重复步骤6～9，直至找到经过缺陷区域的全部射线的最短走时路径。
[0058] 步骤11：将得到的全部最短射线路径作为已知量代入步骤5所述的联合迭代重建算法(SIRT)公式，再次对网格慢度进行求解，得到所有射线在各个小网格内的慢度，然后根据慢度曲线在缺陷区域的突变情况即可判断缺陷大小及分布。如此，如果成像精度不高，可重复步骤2～11，获得更高的成像精度。
[0059] 可变例1
[0060] 步骤1：待测铝板厚度为3mm，成像区域的面积为边长为640mm的正方形区域；标准人造腐蚀坑缺陷的直径为30mm、深为2mm；缺陷中心点位于(320mm， 460mm)。成像区域被划分为128×128个小网格(即128×128个像素)。在成像区域的一侧等间距的设置14个发射EMAT位置，对侧相同位置处设置相同数量的接收EMAT位置，所述的EMAT直径为30mm，相邻EMAT中心间距为45.7mm。
[0061] 步骤2：射频功率放大器选择AG1024，用于激励EMAT，令发射端中第一个发射EMAT全向激发A0模态Lamb波，所有14个接收EMAT依次进行接收，激发频率为290kHz。重复操作，令发射端其他13个发射EMAT均激发过一次A0模态Lamb波。所述的Lamb波为单一的A0模态Lamb波，不含A0以外其他模态Lamb波。
[0062] 步骤3：利用伪Wigner-Ville分布(PWVD)对M×M个EMAT检测波形进行时频分析和模态识别；
[0063] 步骤4：提取所有196个检测波形的分析结果并记录A0模态Lamb波的走时，记为T1～T196。
[0064] 步骤5：使用联合迭代重建算法(SIRT)确定每个网格的慢度(速度的倒数)及缺陷分布：
[0065]
[0066] 其中，Sj为待求的第j个网格的慢度；
[0067] Lij为第i条投影射线在第j个网格中的长度；
[0068] Ti为第i条投影射线的实测走时；
[0069] 可知，n＝128×128＝16384，m＝14×14＝196。
[0070] 计算出的缺陷直径值约为40mm。
[0071] 步骤6：在成像区域建立一个三维直角坐标系，然后进行子网格区域划分，每个子网格区域包含4×4＝16个像素；建立一个30度的初始发射角，利用外推 法逐步对射线路径进行追踪。在外推点区域均匀获取16个像素，得到对应16个cp值，代入二维三次样条插值(2D-CSI)公式，构成16个线性方程，对Lamb波cp的分布情况进行二维曲面拟合，精确求解每一步外推所需的 和
[0072]
[0073] 其中，cp是相速度，x,y表示追踪方向；Ckl(k,l＝0,1,2,3)表示16个待求解系数，16个像素中每一个像素都对应一个x、y值。
[0074] 步骤7：将步骤6得到的cp、 和 代入RT外推公式，完成一步外推：
[0075]
[0076] 其中，yk表示经过k步外推获得的射线i上的离散点沿y轴方向的坐标值；
[0077] yk+1、yk-1分别表示上述离散点前、后两点沿y轴方向的坐标值；
[0078] Δx表示射线外推的固定步长，具体有：
[0079]
[0080] 步骤8：判断新找到的外推点坐标值是否到达或超出子网格区域边界坐标值，若不是，则返回步骤5；若是，则结束对射线i的RT，记录下终点坐标(xf，yf)。
[0081] 步骤9：判断终点坐标(xf，yf)是否到达或接近目标接受位置(即第i对EMAT中的接收EMAT)，若是，则得到射线i的预期最短走时路径；若不是，则每次增加或减少3度，从而改变初始发射角的值，重复步骤6～8，直到找到连接收发EMAT的射线i的最短走时路径。
[0082] 步骤10：重复步骤6～9，直至找到经过缺陷区域的全部射线的最短走时路径。
[0083] 步骤11：将得到的全部最短射线路径作为已知量代入步骤5所述的联合迭代重建算法(SIRT)公式，再次对网格慢度进行求解，得到所有射线在各个小网格内的慢度，然后根据慢度曲线在缺陷区域的突变情况即可判断缺陷大小及分布。
[0084] 如此，重复步骤2～11，以获得更高的成像精度；通常的，缺陷长，宽，深的量化误差小于10％时认为成像精度满足要求。本实施例最后得到的缺陷直径值约为32mm，与实际值吻合良好，量化误差小于10％。
[0085] 以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉
本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应该以
权利要求的保护范围为准。
