CN102636568A - 一种检测混凝土内部缺陷的有限元超声成像方法 - Google Patents

Abstract

本发明公开了一种检测混凝土内部缺陷的有限元超声成像方法，包括将超声发射器和接收器放置在混凝土表面，发射和接收超声波，对采集到的超声信号进行滤波预处理后，用有限元方法对混凝土内部物理参数进行重建，具体步骤为：选择超声传播的模型，信号前处理，有限元迭代重建，以及图像后处理。然后根据图像做结果分析，依据混凝土的结构设计，从而确定混凝土内缺陷的位置、大小和性质。

Description

一种检测混凝土内部缺陷的有限元超声成像方法

技术领域

本发明涉及一种超声成像方法，特别涉及一种检测混凝土内部缺陷的有限元超声成像方法。

背景技术

混凝土由水泥、砂子和石子等基本原材料组成，其经常出现的质量问题包括内部存在不密实或孔洞，外部形成蜂窝麻面、裂缝和损伤层等缺陷，这会严重影响混凝土器件的承载能力和耐久性，甚至造成重大安全事故。超声波检测方法是混凝土无损检测中的一种常规方法，与其它无损检测方法相比，具有指向性好、被测对象范围广、能够检测的深度大、缺陷定位准确、适应性强、检测灵敏度高、成本低、设备体积小、速度快、对人体无害及便于现场使用等优点，因此超声无损检测技术是国内外应用最广泛、使用频度最高且发展较快的一种无损检测技术。

目前，混凝土超声波无损检测方法包括两大类型。一种是直接利用声波在的混凝土中传播的速度、接收波的幅度和频率来判断混凝土的内部缺陷。这仅仅是一种定性的分析技术，若想判断内部缺陷的性质、范围，以及推测强度变化区域就无能为力了。另一种是对接收器接收到的超声信号进行成像，可以提供直观和大量的信息，直接反映探测介质的声学和力学性质，从而直接判断混凝土内部缺陷的存在和属性。超声波成像方法已被应用于混凝土内部缺陷的检测，主要有超声CT法和超声波合成孔径聚焦成像(SAFT)。超声CT对声波走时或衰减进行层析成像，通常假设声波是沿直线传播的，这与实际声波在具有很强散射性的混凝土中的传播有较大出入的。合成孔径聚焦成像是对B扫描或C扫描信号进行合成处理，实现用小孔径传感器来模仿大孔径传感器，获得二维或三维的成像图，其在图像重建过程中假设声速是恒定的，只能重建出结构图，而且当混凝土内部结构较复杂时，重建的结构图不一定准确。

发明内容

本发明解决的技术问题是：提供一种采用有限元方法对混凝土进行超声成像的方法，可以检测混凝土内部的缺陷，克服现有混凝土超声检测方法只能根据超声信号的特征做定性分析的缺点，或进行超声成像时采用超声沿直线传播或声速恒定的假设而成像不准确的缺点。

本发明的技术方案是，一种检测混凝土内部缺陷的有限元超声成像方法，包括以下步骤：

步骤一：将将超声检测仪的多个发射器和发射器放置在混凝土表面，，然后按顺序每次让一个发射器发射超声信号，多个接收器同时接收这个发射器所发射的超声信号，依次进行多次操作，记录不同位置的发射器所发射的超声信号以用于混凝土图像重建；

步骤二：对采集到的超声信号进行滤波处理；

步骤三：对超声信号进行有限元重建，包括选择超声传播的模型，信号前处理，有限元迭代重建以及图像后处理；

步骤四：根据重建后的参数图像，确定混凝土内各位置成分，从而找到缺陷。

所述的一种检测混凝土内部缺陷的有限元超声成像方法，所述的步骤三中，有限元重建方法包括以下过程：

通过采用超声传播模型将所要重建的多组参数Q₁、Q₂等和超声声压P联系起来，即

f(P，Q₁，Q₂，…Q_n)＝0

然后将重建区域分为多个小面元，重建区域内的任一点值由其所在面元顶点的值通过插值确定，并对模型方程进行离散化处理得到方程组，再经重建算法进行求解来求得方程组中的未知数Q₁、Q₂等多组待重建参数，从而得到整个重建区域的图像。

所述的一种检测混凝土内部缺陷的有限元超声成像方法，在有限元求解时，所用的频域超声传播方程为

${\▿}^{2}P\left(r\right)+{k_0}^2(1+O)P\left(r\right)=0$

其中P(r)表示超声声压，k₀为参考介质中的超声的波数，k₀＝ω/c₀，ω为超声的角频率，c₀为在参考介质中的超声传播速度，O为表示物体内需要重建的声学参数的变量

$O=\frac{c_o^2}{c^2}-1+i\frac{{\mathrm{Ac}}_0}{c^2}$

其中c为超声传播速度，A为超声衰减系数，上式可以简写为：

O＝O_R+iO_I

O_R和O_I即为要重建的参数，从而可求得超声传播速度c以及超声衰减系数A。

所述的一种检测混凝土内部缺陷的有限元超声成像方法，在有限元逆向迭代求解时，所用的主要方程为

$(J^{T}J+\mathrm{\λI})\mathrm{\Δ\χ}=J^T(P_s^o-P_s^c)$

其中J为由超声声压对所重建变量的偏导数组成的雅克比矩阵，J^T为J的转置矩阵，λ为重建常数因子，I为单位矩阵，Δχ为迭代中所重建变量的增量，

为接收器接收的声压，

为根据假设的初始值计算得到的声压。

本发明的技术效果在于，首次将混凝土内部的超声声压P(r，t)或P(r，k)和需要重建的n种混凝土内部参数Q₁(r)，Q₂(r)...Q_n(r)用一个超声传播方程联系起来，将混凝土超声成像变成一种基于模型的图像重建方法，即

f(P，Q₁，Q₂，…Q_n)＝0

然后运用有限元方法，首次得到混凝土内部具有物理意义的重建参数分布图。同时采用频域有限元的方法，首次用混凝土内部超声传播速度分布图和超声衰减系数分布图分析混凝土内部的结构，从而得到缺陷或者目标的位置。其突出的优点在于：

1.相比于普通的、根据超声信号进行定性分析的方法，它直接提供混凝土内部重建参数的分布图，因此更为直观，简单。

2.它消除了普通超声混凝土成像方法中超声按直线路径传播或传播速度处处相同的假设，因此在判断缺陷空间结构和位置上也更为准确。

3.相对于其它超声成像方法，本方法适用范围广，还可以用在除混凝土之外的其它固体内部结构探测中。

4.相对于普通超声成像方法，本方法对超声发射器和接收器位置没有太多要求，使用灵活；而且通过采用不同超声传播模型，可以重建不同的声学参数，功能强大。

附图说明

图1为本发明的步骤流程图；

图2为混凝土试件结构示意图以及超声发射器和接收器分布示意图；

图3为有限元方法中对成像区域的网格划分图；

图4为有限元重建流程图

图5为频域有限元算法流程图

图6为信号预处理后的超声波形图

图7为图6所对应的傅里叶变换后的频谱图；

图8为重建后的超声传播速度分布图；

图9为重建后的超声衰减系数分布图。

具体实施方式

下面结合具体实例，对本发明技术方案做进一步说明。在本实例中，我们要对一个含有金属波纹管的预应力混凝土梁试件进行横截面成像，得到超声传播速度分布图和超声衰减系数分布图，从而判断混凝土内部的结构和成分。叙述顺序如图1发明流程图所示。

步骤100：发射和接收超声波

首先，假设在混凝土表面有S个超声发射器位置，N个超声接收器位置。首先第一个位置的超声发射器发出超声，所述N个位置上的超声接收器接收一系列超声信号，记为P₁(t)、P₂(t)...P_N(t)；然后第二个位置的超声发射器发出超声，所述N个位置的超声接收器再次接收一系列超声信号，记为P_N+1(t)、P_N+2(t)...P_2N(t)；按照这种方式，总共收集S*N个接收器位置的信号。

图2为本例中混凝土试件的结构图，其长为600mm，高为504mm。试件的上下两侧分别设有18个超声发射器位置，18个超声接收器位置。按照上述方法收集18*18＝324个超声信号，信号的采样频率为2MHz，采样长度为2000sp。

步骤200：信号预处理

根据需要，对接收到的所有超声信号进行去噪、截取等滤波处理，以便于下步计算。

步骤300：有限元重建

如图3所示，它包括模型和重建算法的选择，信号前处理，有限元迭代重建，以及图像后处理四部分。这一步骤是本方法的重点。

(1)超声传播模型的选择

首先有限元计算是一种基于模型的数值计算方法，因此首先要选择所基于的模型，同时也就确定了需要重建的参数。在有限元超声成像中，就是要找到一个超声传播方程，将所要重建的多组参数Q₁、Q₂等和超声声压P联系起来，即

f(P，Q₁，Q₂，…Q_n)＝0               (1)

在本例中，模型采用频域下超声的传播模型

${\▿}^{2}P\left(r\right)+{k_0}^2(1+O)P\left(r\right)=0---\left(2\right)$

其中波数k₀＝ω/c₀，ω为超声的角频率，c₀为在参考介质中的超声传播速度，O为表示混凝土内部性质的参数

$O=\frac{c_o^2}{c^2}-1+\frac{i}{k_0}\frac{{\mathrm{Ac}}_0}{c^2}---\left(3\right)$

其中c为超声传播速度，A为超声衰减系数。我们可以简写为：

$O=O_R+i\frac{O_I}{k_0}---\left(4\right)$

因为超声波P可以分为入射波P_inc和散射波P_s，P(r)＝P_inc(r)+P_s(r)，而入射波满足Helmholtz方程，即

所以模型可变为：

${\▿}^2{P}_s+k_0^2(1+O)P_s=-k_0^2{\mathrm{OP}}_{\mathrm{inc}}---\left(5\right)$

(2)频域有限元算法

有限元主要思想为将重建区域分为很多小面元，如图4。重建区域内的任一点值由其所在面元顶点的值通过插值确定，然后用数学手段就可以将所用的模型方程变为一组方程组，即进行离散化处理。这样，通过反演算法，就可以求得该方程组中的未知数，从而得到整个重建区域的图像。

在本例中采用频域有限元算法，采用式(5)作为超声传播模型。图3为对整个成像区域进行网格划分后的结果，其中超声发射器和接收器的位置如图中所示。

本频域有限元算法表述如下：

首先将成像区域中任一一点的散射波声压表示成如下式子：

$P_s(r,k)=\underset{i=1}{\overset{3}{\mathrm{\Σ}}}{P}_{s,i}\left(k\right){\mathrm{\ψ}}_i\left(r\right)---\left(6\right)$

其中P_s，i(k)为其所在面元的上第i个节点在频率为k时的散射波声压值，ψ_i(r)为第i个节点的插值函数，这个插值函数只和网的划分有关。上面i＝1，2，3是因为采用三角形面元划分，每个面元有3个节点。同理，方程内的O_R、O_I、P_inc都需要用这种方法插值表示。需要注意的是，可以对成像区域用不同的网进行划分，让不同的参数定义在不同的网上，在这里我们只取简单的做法。

在插值完毕后，通过将(5)式用(6)式进行代换，然后方程两边同时乘以某一个节点j的插值函数，再在全成像区域上对方程两边进行积分，(5)式就可以表示成如下的形式：

其中，在进行全成像区域积分的时候，用到了格林函数积分公式，引入了如下的边界条件：

其中， $\mathrm{\α}=\frac{-\mathrm{ik}-3/2\mathrm{\ρ}+i3/8{\mathrm{k\ρ}}^2}{1-i/\mathrm{k\ρ}};$ $\mathrm{\β}=\frac{-i/2k{\mathrm{\ρ}}^2}{1-i/\mathrm{k\ρ}},$ N为整个网的节点数，而入射波P_inc可以通过格林函数 $G\left({r,r}_0\right)=\∫\frac{e^{\mathrm{ik}|r-r_0|}}{\sqrt{|r-r_0|}}\mathrm{ds}$ 求出来：

P_inc＝∫F(r₀)G(r，r₀)ds     (9)

其中r₀代表发射器的位置，为已知量。积分过后，(7)式就可以表示成一个方程组：

[A]{P_s}＝[B]{P_inc}          (10)

其中

$B_{\mathrm{ij}}=k_0^2<\underset{k}{\mathrm{\&Sigma;}}{O}_{R,k}{\mathrm{\&psi;}}_k{\mathrm{\&psi;}}_j{\mathrm{\&psi;}}_i>+{\mathrm{ik}}_0<\underset{l}{\mathrm{\&Sigma;}}{O}_{I,l}{\mathrm{\&psi;}}_l{\mathrm{\&psi;}}_j{\mathrm{\&psi;}}_i>$

{P_s}＝{P_s，1，P_s，2，…，P_s，N}^T

{P_inc}＝{P_inc，1，P_inc，2，…，P_inc，N}^T

在(10)式中，P_inc、[B]已知，需要重建的参数O_R和O_I在[A]中，P_s未知。因为已知接收器处的真实声压值

它是成像区域内的重建参数真实分布

和

的函数。为了重建[A]中的未知参数，可以将接收器处的实测散射波写成如下的形式：

$P_s({\tilde{O}}_R,{\tilde{O}}_I)=P_s(O_R,O_I)+\frac{\&PartialD;P_s}{\&PartialD;O_R}\mathrm{\&Delta;}{O}_R+\frac{\&PartialD;P_s}{\&PartialD;O_I}\mathrm{\&Delta;}{O}_I+\&CenterDot;\&CenterDot;\&CenterDot;---\left(11\right)$

其中P_s(O_R，O_I)为在参数分布为O_R、O_I时，根据(10)式计算得到的值；另外 $\mathrm{\&Delta;}{O}_R={\tilde{O}}_R-O_R,$ $\mathrm{\&Delta;}{O}_I={\tilde{O}}_I-O_I.$

考虑到接收器的数据有M组，其中包括所有的发射器位置S，接收器位置ND，以及频率数目NF，即M＝S*ND*NF；而需要重建的参数O_R和O_I各有K和L，那么在一阶近似下上式可以改写为：

$\mathrm{J\&Delta;\&chi;}=P_s^o-P_s^c---\left(12\right)$

其中J为雅克比矩阵，

$\mathrm{\&Delta;\&chi;}=\left\{\begin{array}{c}\mathrm{\&Delta;}{O}_{R,1}\\ \mathrm{\&Delta;}{O}_{R,2}\\ .\\ .\\ .\\ \mathrm{\&Delta;}{O}_{R,K}\\ \mathrm{\&Delta;}{O}_{I,1}\\ \mathrm{\&Delta;}{O}_{I,2}\\ .\\ .\\ .\\ \mathrm{\&Delta;}{O}_{I,L}\end{array}\right\},$ $P_s^o=\left\{\begin{array}{c}{P}_{s,1}^o\\ P_{s,2}^o\\ .\\ .\\ .\\ P_{s,M}^o\end{array}\right\},$ $P_s^c=\left\{\begin{array}{c}{P}_{s,1}^c\\ P_{s,2}^c\\ .\\ .\\ .\\ P_{s,M}^c\end{array}\right\}$

由于采用了一阶近似，这里把

改写为P^o，因为在近似情况下它并不是测量到的声压值，而只是一种接近的情况；而P_s(O_R，O_I)为在参数分布为O_R、O_I时，根据(10)式计算得到的值，因此标记为P^c。散射波对要重建参数O_R的偏导

可以根据(10)式求导得到

$\left[A\right]\left\{\frac{{\&PartialD;P}_s}{\&PartialD;O_R}\right\}=\left[\frac{\&PartialD;B}{\&PartialD;O_R}\right]\left\{P_{\mathrm{inc}}\right\}-\left[\frac{\&PartialD;A}{\&PartialD;O_R}\right]\left\{P_s\right\}---\left(11\right)$

同理，对于也有

$\left[A\right]\left\{\frac{{\&PartialD;P}_s}{\&PartialD;O_I}\right\}=\left[\frac{\&PartialD;B}{\&PartialD;O_I}\right]\left\{P_{\mathrm{inc}}\right\}-\left[\frac{\&PartialD;A}{\&PartialD;O_I}\right]\left\{P_s\right\}---\left(12\right)$

然后对(12)式作变换，两边同时乘以雅克比矩阵的转置J^T，并且引入一个重建常数因子λ，得到

$(J^{T}J+\mathrm{\&lambda;I})\mathrm{\&Delta;\&chi;}=J^T(P_s^o-P_s^c)---\left(13\right)$

其中I为单位矩阵。

求解(13)式，就可以得到ΔO_R和ΔO_I，从而解得和

的图像。由于

和超声传播速度有关，而

和超声衰减系数以及超声传播速度有关，因此由和

的值就可以得到超声传播速度分布图以及超声衰减系数分布图。

然而这只是一步计算的结果，它不一定是重建参数的真实分布，如果图像效果不是很好，为了得到更好的结果，可以把这组值当做下一步的初始分布再代入(10)式中进行计算得到新的P_s(O_R，O_I)，并且用(11)式和(12)式得到雅克比矩阵J，然后再都代入到(13)式计算新的

和

这样，一直用(10)式和(13)式迭代，直到得到满意的图像为止。其算法结构流程图如图5。

(3)信号前处理

下面继续以图2所示试件为例，说明本方法的应用。

信号前处理要根据所用的有限元算法进行。因为本例中用的是频域有限元算法，因此先要对采集到的信号作傅里叶变换。图6为图2中超声发射器在S5位置且超声接收器在D12位置接收到的超声信号做了信号预处理之后的波形图。其中波形图共有2000个采集点，时间长度为0.5毫秒。从图中可以看出，在0.13毫秒处有一个非常明显的峰，这个峰正是超声接收器接收到的超声发射器发射出的声脉冲传过来的信号。而在这个信号之后，有一系列的信号的影响，这是波纹管在声脉冲作用下产生出来的波纹管内部的信号。

图7为对图6做了傅里叶变换后的频谱图。图中共有1000个频率，最高频率达到1MHz其中从图中可以看出，所采集的信号频率集中在200kHz一下。因此我们从这1000个频率中，等间隔的在0到200kHz中挑选10个频率，用来进行下一步的有限元重建。

即在本方法中，信号前处理即对处理后的采集信号作傅里叶变换，然后挑选其中的某些频率信号用于下一步的计算。

(4)有限元迭代重建

信号前处理后，就可以用本方法中的频域有限元方法进行迭代重建了。本例中，共有18个超声发射器位置，18个超声接收器位置，信号前处理中选择了10个频率，因此共有已知数数量为18*18*10＝3240个。用来划分成像区域的网共有节点1924个，如图4所示。因为超声传播速度和超声衰减系数都是定义在这个网上，因此共有未知数的数量1924*2＝3848个。

整个重建过程共用了3次迭代，用时大约1个小时。

(5)图像后处理

得到了有限元重建结果，即所用网格的节点上的重建参数值后，就可以查看图像了，并对图像做一些简单滤波处理。本例中处理后的超声传播速度分布图和超声衰减系数分布图如图8和图9所示。

步骤400：结果分析

通常，混凝土中缺陷和其周围混凝土的性质(力学参数或声学参数)存在较大的差异，本发明重建得到的参数分布图可以用于检测缺陷的存在，并做定量分析。

如图8和图9所示，无论是在超声传播速度分布图还是在超声衰减系数分布图中，都可以明显地看到混凝土试件中金属波纹管的位置及其内部结构。从两幅图中量出金属波纹管的直径大小分别为55mm和52mm，和实际的直径50mm非常接近，因此可以说本方法很好的得到了混凝土中的结构。更为重要的是，从超声速度传播分布图中可以看到，背景混凝土的超声传播速度平均为4000m/s，金属波纹管内超声的速度分为两部分，其中靠波纹管左侧、平均值约为4800m/s部分可以认为是钢绞束，而另一侧平均值约为3400m/s部分可以认为是波纹管压浆不密实区。而在超声衰减系数分布图中，可以看到背景的衰减系数为0.15，而波纹管内的衰减系数约为0.03(右侧的衰减系数比左侧略大)，这是因为波纹管中钢绞束对声波几乎不衰减，而不密实区的衰减系数比背景混凝土的衰减系数小。

通过对所述示例混凝土试件进行破损检测，上述图像重建的结论得到了验证。

Claims (4)

1.一种检测混凝土内部缺陷的有限元超声成像方法，其特征在于，包括以下步骤：

步骤一：将超声检测仪的多个发射器和发射器放置在混凝土表面，，然后按顺序每次让一个发射器发射超声信号，多个接收器同时接收这个发射器所发射的超声信号，依次进行多次操作，记录不同位置的发射器所发射的超声信号以用于混凝土图像重建；

步骤二：对采集到的超声信号进行滤波处理；

步骤三：对超声信号进行有限元重建，包括选择超声传播的模型，信号前处理，有限元迭代重建以及图像后处理；

步骤四：根据重建后的参数图像，确定混凝土内各位置成分，从而找到缺陷。

2.根据权利要求1所述的一种检测混凝土内部缺陷的有限元超声成像方法，其特征在于：所述的步骤三中，有限元重建方法包括以下过程：

通过采用超声传播模型将所要重建的多组参数Q₁、Q₂等和超声声压P联系起来，即

f(P，Q₁，Q₂，…Q_n)＝0

然后将重建区域分为多个小面元，重建区域内的任一点值由其所在面元顶点的值通过插值确定，并对模型方程进行离散化处理得到方程组，再经重建算法进行求解来求得方程组中的未知数Q₁、Q₂…Q_n多组待重建参数，从而得到整个重建区域的图像。

3.根据权利要求2所述的一种检测混凝土内部缺陷的有限元超声成像方法，其特征在于：在有限元求解时，所用的频域超声传播方程为

${\&dtri;}^{2}P\left(r\right)+{k_0}^2(1+O)P\left(r\right)=0$

其中P(r)表示超声声压，k₀为参考介质中的超声的波数，k₀＝ω/c₀，ω为超声的角频率，c₀为在参考介质中的超声传播速度，O为表示物体内需要重建的声学参数的变量

$O=\frac{c_o^2}{c^2}-1+i\frac{{\mathrm{Ac}}_0}{c^2}$

其中c为超声传播速度，A为超声衰减系数，上式可以简写为：

O＝O_R+iO_I

O_R和O_I即为要重建的参数，从而可求得超声传播速度c以及超声衰减系数A。

4.根据权利要求2所述的一种检测混凝土内部缺陷的有限元超声成像方法，其特征在于：在有限元逆向迭代求解时，所用的方程为

$(J^{T}J+\mathrm{\&lambda;I})\mathrm{\&Delta;\&chi;}=J^T(P_s^o-P_s^c)$

其中J为由超声声压对所重建变量的偏导数组成的雅克比矩阵，J^T为雅可比矩阵J的转置矩阵，λ为重建常数因子，I为单位矩阵，Δχ为迭代中所重建变量的增量，

为接收器接收的声压，为根据假设的初始值计算得到的声压。

Images