CN113252779B - 一种提高缺陷反演成像质量的方法 - Google Patents

Abstract

本发明提供一种提高缺陷反演成像质量的方法，首先对无缺陷板中的导波信号进行傅里叶变换取得幅值。再对有缺陷的板中导波信号进行傅里叶变换取得幅值，将两者幅值相减，实现抵消掉无损铝板边界及原生源影响的效果，从而保留下由缺陷产生的散射波幅值。然后，利用矩量法将算子方程离散化为矩阵方程，并构建超定方程来对矩阵方程进行求解。求解过程中利用截断最小二乘法实现矩阵方程的正则化，通过改进的奇异值曲线方法选取正则化参数，克服广义交叉验证方法中无法准确选择正则化参数的问题。最后，通过截去部分奇异值来减少信号噪声，降低噪声干扰，从而有效提高缺陷反演的准确性。

Description

一种提高缺陷反演成像质量的方法

技术领域

本发明涉及无损检测领域中的超声导波检测技术，具体涉及一种提高缺陷反演成像质量的方法。

背景技术

传统超声波检测方法有效检测区域仅仅为传感器下方声束扫描的有限区域，对于探头难以到达的位置无法实现检测，与之相对，超声导波在板中仅需由单点激励，在呈环形阵列分布的多个接收换能器上进行接收，便可实现对整个板材缺陷的快速扫查。

目前来说对于板材成像，大部分都是用假设导波在板中沿直线传播的方法，忽略了缺陷处发生的散射、衍射及模态转换等物理现象；而针对介质密度变化不大的板材提出的射线追踪的方法，尽管是引进了射线追踪法来追踪声波传播路径，再根据路径对重建算法进行了改进，但是其忽略了衍射，同样当出现衍射这一物理现象时会导致失效。所以射线追踪法对于导波沿直线传播的改进仍不具备普适性。

在计算成像解时由于构建的方程为超定方程且采集的数据往往伴随噪声污染，所以正则化对获得精确解尤为重要，目前主流、处于核心地位的对处于噪声干扰下的系统求解的正则化方法是Tikhonov正则化方法，但是此方法仅对病态矩阵修正，某些情况下无法达到理想去噪效果。

发明内容

针对上述内容已有问题，本发明提供一种通过导波幅值来计算的基于矩量法的提高缺陷检测成像质量的方法。

实现本发明目的的技术方案是：一种提高缺陷反演成像质量的方法，对换能器接收到的有损铝板的导波与无损铝板的导波分别进行傅里叶变换，取得傅里叶变换后产生的频域图中中心频率对应的幅值，并且将有缺陷铝板与无缺陷铝板所对应的不同幅值进行相减，获得的差值当作缺陷作为次生源发射的散射波幅值，用矩量法构建超定方程来表示各离散点作为次生源发射lamb波的叠加作用，最终对此成像超定方程实现求解，在本次计算中将铝板划分为1793个单元点，每个单元点之间距离为2mm，由于要计算每个单元点上的未知数，所以构建的超定方程组，方程数需大于1793个。

本发明采用矩量法构建表示铝板中各离散点介质声学特性的参数的方程组，通过picard准则对方程进行正定性判断，将用于正则化的截断奇异值法与最小二乘法结合实现方程求解，确定各点的声速(折射系数)特性参数，进而实现成像。

其中，具体装置包括：有缺陷铝板、无缺陷铝板、高能超声测试系统、超声导波换能器、信号发生器、示波器；在所述两种铝板上放置换能器。

一种提高缺陷反演成像质量的方法，包括：

S1：获取换能器在无损铝板上接收到的无损导波信号，及换能器在具有缺陷的有损铝板上接收到的有损导波信号；

S2：无损导波信号和有损导波信号分别进行傅里叶变换，对获取的幅值进行相减，得到缺陷作为次生源发射的散射波幅值，构成散射矩阵p^(s)；

S3：采用矩量法构建缺陷中各离散点作为次生源发射lamb波叠加作用的超定方程；

S4：通过picard准则对超定方程进行正定性判断，将用于正则化的截断奇异值法与最小二乘法结合实现超定方程求解，获得缺陷中各离散点的介质声学特性参数，实现缺陷成像。

求解过程中的截断系数由L曲线法获得。

进一步的，所述步骤S3的具体步骤包括：

S31：建立超声全场与散射场算子方程L(f)＝g，其中L为线性算子，f为用来求解的未知函数，g为已知的缺陷次生源函数；通过波动方程建立超声全场与散射场方程L(f)＝g，由算子方程来表示声场声压分布。

S32：将算子方程离散化为矩阵方程，离散点应足够多，使此方程为超定方程，通过Picard准则对其进行受噪声污染程度判断。

此时，在算子方程L(f)＝g中，g函数为散射矩阵p^(s)，f函数为用来求解的表示物体内部各离散点的介质声学特性的参数矩阵O。

进一步的，所述步骤S31的具体步骤包括：

由运动方程、连续方程、物态方程推出适用于铝板的非齐次波动方程：

其中，为/>点的折射系数，/>c₀为在均匀介质中声波的传播速度；/>为非均匀介质内的折射系数/>变化的函数，即/>n_δ表示折射系数的变化，k(r)为波数，/>在非均匀介质中/> 为原均匀介质的波数，/>为声压；

非均匀介质中超声波作用产生的压力场表示为超声波穿过均匀介质产生的入射场/>和穿过非均匀介质所产生的散射场/>的和，用/>来表示；

入射场表示为：

散射场表示为：

添加关于距离的格林函数，得散射场算子方程与全场算子方程：

优选的，所述步骤S31的具体步骤包括：

由运动方程、连续方程、物态方程推出适用于铝板的二维声波方程：

此处用拉普拉斯算符代替/>在非均匀介质中/>根据得到非均匀介质中/>点的折射系数/>

c₀为在均匀介质中声波的传播速度，把看作非均匀介质内的折射系数/>变化的函数，即/>n_δ表示折射系数的变化。若同时考虑声速特性参数与衰减特性参数，则未知函数可表示为：

由于本成像方法中仅考虑导波的声速特性参数，不考虑衰减特性参数，所以设得非齐次波动方程：

进一步有全场方程及散射场方程：、

进一步的，步骤S32的具体步骤包括：

把要探究的区域划分成的网格，本方法中设n＝1793，把狄拉克函数δ作为矩量法的检验函数，

矩量法中用离散和来代替连续积分，存在缺陷发生散射，在此路径的和将会较大，但若该点没有散射体，求和便会小很多。二维情况下，格林函数能够被第一类汉克尔函数表示：/>因为此处是对小单元进行的分析，所以在小单元内做内接圆，进而获得

m＝1,2,3,…M,表示接收信号的探头。为了简化表示，将矩阵和向量引入：

C＝c_pq为N×N维的矩阵，是未知函数构成的矩阵，D＝d_mq是M×N维矩阵，得到矩阵方程组：

p^(t)＝p⁽ⁱⁿ⁾+COp^(t)

p^(s)＝DOp^(t)

进一步的，步骤S4的具体步骤包括：

假设超定方程为Ax＝b，即或/>

第一步，建立增广矩阵[A,b]，然后对增广矩阵进行奇异值分解；

第二步，选取截断参数k，

k≤min(n,rank(A,b)),σ_k＞σ_k+1且V₂₂≡(v_n+1,k+1,…v_n+1,n+1)≠0；

第三步，记q＝n-k+1，对第一步中奇异值分解获得的V进行分块，划分为四部分：

最后通过第二步获取的截断参数，实现对给定数据拟合质量与范数最小的权衡；

通过完全最小二乘法获得的解的范数可表示为：

体现拟合程度的残差范数可表示为：

完全最小二乘的解为

进一步的，所述步骤S4的具体步骤包括：

获得缺陷中各离散点的介质声学特性的参数的值，由此值作为构建的三维图的纵坐标z，离散后各点的行列位置对应x、y坐标，在matlab中根据离散点坐标以及colormap里面的颜色表自动确定颜色，进而实现缺陷成像。

本发明基于矩量法的超声导波信号差值成像方法，具有以下有益效果：

相较于目前普遍使用的仅仅用到衍射波或反射波的方法，本方法由于实现了对衍射波、折射波、反射波的利用，因此能够得到更好的成像效果，且能够适用于不同深度板材缺陷，不像衍射波或者折射波成像方法对板材缺陷深度有要求限制。最终通过计算获得的声速(折射系数)特性参数实现成像；

对适用于Tikhonov正则化参数选取的奇异值曲线方法进行改造，使其适用于截断最小二乘法截断参数的选取，克服了广义交叉验证方法中存在的无法准确选择哪个正则化参数的问题。且不需要像广义交叉验证方法那样对数据噪音情况有所了解或者估计；

本发明将矩量法应用到超声导波领域，基于换能器测得的数据，将波动方程、矩量法、截断最小二乘法等进行了详细推导，对成像效果进行了数值模拟，且成像位置准确、能体现内部构造。

附图说明

图1为发明中待检测铝板缺陷位置图。

图2为发明中用于选取截断参数位置的奇异值曲线图。

图3为发明中用于判断模型受污染程度的Picard准则图。

图4为发明中对仿真数据成像图。

具体实施方式

为了使本发明的技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明具体实施方式进行描述。应当理解，此处描述的具体实施例仅用以解释本发明，并不用于限定本发明。

本发明公开的一种提高缺陷反演成像质量的方法。通过矩量法，将通过波动方程构建的全场、散射场方程组由算子方程转换为代数方程。

求解过程中引进截断完全最小二乘法来选择正则化参数，进而解决求解震荡发散的问题。在代数方程不适定性较弱时，截断完全最小二乘法方法对Tikhonov正则化的方法优势并不明显，但是当噪音的影响较大，方程具强不适定性时，本方法中选择的截断完全最小二乘法相比较Tikhonov正则化方法具备明显的优势，能够达到提高成像的分辨率的目的。

本发明将换能器分别在无缺陷板与有缺陷板中接收到的信号进行傅里叶变换，通过变换后幅值相减的方式来保留仅有缺陷作为次生源发出的信号。然后通过矩量法将算子方程转化为矩阵方程，构建超定方程求解，求解过程通过截断最小二乘法实现方程的正则化，完全截去部分奇异值实现噪声污染的减弱。

本发明通过矩量法实现对所建立的超定方程的求解，借助求解出的折射系数实现缺陷成像检测。

根据本申请的一个实施例，建立仿真，参考图1。在距铝板中心水平距离29mm处，存在一半径为10mm的缺陷，缺陷深度2mm，以铝板中心为圆心，等距的设置64个用来放置探头的位置，使其呈环形阵列分布，环形阵列半径50mm。lamb波激励频率为250khz，无损板中波速2324.7m/s。

S1：对在无损铝板上放置换能器接收到的波及在有损铝板上放置换能器接收到的波分别进行傅里叶变换，进而获得中心频率对应的振幅，并求两者差值。由于选择了64个位置点分别放置探头激励，每次激励对应63个不同位置点分别放置探头接收，所以获得4032个振幅差值，进而构建4032×1的散射矩阵p^(s)。

S2:通过波动方程建立超声全场与散射场方程L(f)＝g，由算子方程来表示声场声压分布。

S3：将算子方程经过离散化，成为矩阵方程。离散点应足够多，使此方程为超定方程，本超定方程组中具有1793个未知数，方程组中有4032个方程，通过Picard准则对其进行受噪声污染程度判断，方程组Ax＝b的傅里叶系数<u_i，b>趋于零的速度在平均意义下快于矩阵A的奇异值σ_i趋于零的速度，Picard准则判断图见于图3。

在计算过程中，在铝板中心以半径100mm做圆，作感兴趣区域。离散点相互之间距离2mm。

S4：借助截断最小二乘法对超定方程进行求解，获得折射系数k＝90。求解过程中的截断系数由奇异值曲线法获得。

以下为对上述步骤进行详细描述

步骤S1：对换能器接收到的无损铝板波及有损铝板的波进行傅里叶变换获得幅值求差值。

分别对有缺陷与无缺陷铝板采取一发，其余63个探头位置分别放置接收探头进行接收的数据采集模式，获得两组数据，每组数据为4032个波形数据，对每组数据进行傅里叶变换，提取中心频率振幅，并且两组数据相减，获取4032个散射波数据幅值，进而构建4032×1的散射矩阵p^(s)。

步骤S2:通过波动方程建立超声全场与散射场方程L(f)＝g，由算子方程来表示声场声压分布。

由运动方程、连续方程、物态方程推出适用于铝板的非齐次波动方程：

为/>点的折射系数，/>c₀为在均匀介质中声波的传播速度，把/>看作非均匀介质内的折射系数/>变化的函数，即/>n_δ表示折射系数的变化，k(r)为波数，/>在非均匀介质中/> 为原均匀介质的波数，/>为声压。

通常非均匀介质中超声波作用产生的压力场可以表示为超声波穿过均匀介质产生的入射场/>和穿过非均匀介质所产生的散射场/>的和，可用/>来表示。据此，入射场/>可表示为：

散射场可表示为：

添加关于距离的函数——格林函数，可得散射场方程与全场方程：

步骤S3：将算子方程经过离散化，成为矩阵方程。离散点应足够多，使此方程为超定方程，通过Picard准则对其进行受噪声污染程度判断。

把要探究的区域划分成的网格，可得：

q＝1,2,3,…n，所以全场方程可表示为

把狄拉克函数δ作为检验函数(权函数)，与上式两边做内积，则区域内任意r_p点可表达为：

p为划分基函数数目，即单元数，p＝1,2,3,…n，根据检验函数的性质，则有：

矩量法中用离散和来代替连续积分，存在缺陷发生散射，在此路径的和将会较大，但若该点没有散射体，求和便会小很多。二维情况下，格林函数能够被第一类汉克尔函数表示：/>因为此处是对小单元进行的分析，所以在小单元内做内接圆，来获得以下式子：

其中a是内接圆半径，R_pq＝|r_p-r_q.|，J₁是第一类贝塞尔函数，是n阶第二类汉克尔函数，J_n是n阶第一类贝塞尔函数，N_n是n阶第一类诺伊曼函数，即第二类贝塞尔函数，简化可得，

则全场可表示为

同样，对与散射场进行离散化，得到下式：

m＝1,2,3,…M,表示接收信号的探头。为了简化表示，将矩阵和向量引入：

C＝c_pq为N×N维的矩阵，是未知函数构成的矩阵，D＝d_mq是M×N维矩阵。

p^(t)＝p⁽ⁱⁿ⁾+COp^(t)

p^(s)＝DOp^(t)

上述矩阵便为构建的超定方程，必要条件是划分单元格数小于4020，4020由激励探头位置数64与每次激励后接收探头的位置数63相乘计算获得。

步骤S4：借助截断最小二乘法对超定方程进行求解，获得折射系数。求解过程中的截断系数由L曲线法获得。

假设上面方程为Ax＝b，即或/>

第一步建立增广矩阵[A,b]，然后对增广矩阵进行奇异值分解。

第二步，选取截断参数k，目的是降低小的奇异值对结果的影响。

k≤min(n,rank(A,b)),σ_k＞σ_k+1且V₂₂≡(v_n+1,k+1,…v_n+1,n+1)≠0

第三步，记q＝n-k+1，对(1)中奇异值分解获得的V进行分块，划分为四部分：

最后对于第二步获取截断参数，由奇异值概念可知，矩阵中隐含的信息，重要性与奇异值大小存在正相关关系，且增广矩阵[A,b]可分解成 的形式。所以我们建立如图2所示的关于增广矩阵/>奇异值的奇异值曲线来对截断参数进行选取。得出完全最小二乘的解为/>成像表示于图4，选取截断参数k＝90。

应当理解的是，对本领域普通技术人员来说，可以根据上述说明加以改进或变换，而这些改进和变换都应在保护范围之内。

Claims (4)

1.一种提高缺陷反演成像质量的方法，其特征在于，包括：

S1：获取换能器在无损铝板上接收到的无损导波信号，及换能器在具有缺陷的有损铝板上接收到的有损导波信号；

S2：无损导波信号和有损导波信号分别进行傅里叶变换，对获取的幅值进行相减，得到缺陷作为次生源发射的散射波幅值，构成散射矩阵p^(s)；

S3：采用矩量法构建缺陷中各离散点作为次生源发射lamb波叠加作用的超定方程；

S4：通过picard准则对超定方程进行正定性判断，将用于正则化的截断奇异值法与最小二乘法结合实现超定方程求解，获得缺陷中各离散点的介质声学特性参数，实现缺陷成像；

所述步骤S3的具体步骤包括：

S31：建立超声全场与散射场算子方程L(f)＝g，其中L为线性算子，f为用来求解的未知函数，g为已知的缺陷次生源函数；

S32：将算子方程离散化为矩阵方程，使矩阵方程为超定方程；

此时，在算子方程L(f)＝g中，g函数为散射矩阵p^(s)，f函数为用来求解的表示物体内部各离散点的介质声学特性的参数矩阵O

所述步骤S31的具体步骤包括：

由运动方程、连续方程、物态方程推出适用于铝板的非齐次波动方程：

其中，为/>点的折射系数，/>c₀为在均匀介质中声波的传播速度；/>为非均匀介质内的折射系数/>变化的函数，即/>n_δ表示折射系数的变化，k(r)为波数，/>在非均匀介质中/> 为原均匀介质的波数，/>为声压；

非均匀介质中超声波作用产生的压力场表示为超声波穿过均匀介质产生的入射场/>和穿过非均匀介质所产生的散射场/>的和，用/>来表示；

入射场表示为：

散射场表示为：

添加关于距离的格林函数，得散射场算子方程与全场算子方程：

其中，所述步骤S31的具体步骤包括：

由运动方程、连续方程、物态方程推出适用于铝板的二维声波方程：

此处用拉普拉斯算符代替/>在非均匀介质中/>根据得到非均匀介质中/>点的折射系数/>

c₀为在均匀介质中声波的传播速度，把看作非均匀介质内的折射系数/>变化的函数，即/>n_δ表示折射系数的变化。若同时考虑声速特性参数与衰减特性参数，则未知函数可表示为：

由于本成像方法中仅考虑导波的声速特性参数，不考虑衰减特性参数，所以设得非齐次波动方程：

进一步有全场方程及散射场方程：

2.根据权利要求1所述的一种提高缺陷反演成像质量的方法，其特征在于，所述步骤S32的具体步骤包括：

划分的网格区域，可得：

q＝1，2，3，…n，全场算子方程表示为：

把狄拉克函数δ作为检验函数，与上式两边做内积，则区域内任意r_p点可表达为：

p为划分基函数数目，即单元数，p＝1，2，3，…n，根据检验函数的性质，则有：

二维情况下，格林函数能够被第一类汉克尔函数表示：

在离散点小单元内做半径为a的内接圆，a的值为1mm，获得以下式子：

其中a为离散小单元的内接圆半径，R_pq＝|r_p-r_q|为表示离散点p与q之间距离的矩阵，点p、q均用来表示检测区域被离散后存在的任意一点，J₁是第一类贝塞尔函数，是n阶第二类汉克尔函数，J_n是n阶第一类贝塞尔函数，N_n是n阶第一类诺伊曼函数，即第二类贝塞尔函数，简化可得：

则全场算子方程表示为：

同样，对于散射场进行离散化，得到下式：

m＝1，2，3，…M表示接收信号的探头代号，为了简化表示，将矩阵和向量引入：

C＝c_pq为N×N维的矩阵，是未知函数构成的矩阵，D＝d_mq是M×N维矩阵；

p^(t)＝p⁽ⁱⁿ⁾+COp^(t)

p^(s)＝DOp^(t)

上述矩阵便为构建的超定方程，必要条件是划分单元格数小于63＊64。

3.根据权利要求2所述的一种提高缺陷反演成像质量的方法，其特征在于，所述步骤S4的具体步骤包括：

假设超定方程为Ax＝b，即或/>

第一步，建立增广矩阵[A，b]，然后对增广矩阵进行奇异值分解；

第二步，选取截断参数k，

k≤min(n,rank(A,b)),σ_k>σ_k+1且V₂₂≡(v_n+1，k+1，…v_n+1，n+1)≠0；

第三步，记q＝n-k+1，对第一步中奇异值分解获得的V进行分块，划分为四部分：

V₁₁∈R^n×k,V₁₂∈R^n×q,v₂₁∈R^l×k,V₂₂∈R^l×q

最后通过第二步获取的截断参数，实现对给定数据拟合质量与范数最小的权衡；

通过完全最小二乘法获得的解的范数可表示为：

体现拟合程度的残差范数可表示为：

完全最小二乘的解为

4.根据权利要求1所述的一种提高缺陷反演成像质量的方法，其特征在于，所述步骤S4的具体步骤包括：

获得缺陷中各离散点的介质声学特性的参数的值，由此值作为构建的三维图的纵坐标z，离散后各点的行列位置对应x、y坐标，在matla b中根据离散点坐标以及colormap里面的颜色表自动确定颜色，进而实现缺陷成像。