CN1908651A - 混凝土超声层析成像算法 - Google Patents

Abstract

本发明公开了一种混凝土超声层析成像算法，提出了塔式ART算法，适用于在超声无损检测现场对被测混凝土进行实时二维反演成像，从而准确反映混凝土内部结构。塔式ART算法将网格逐步分块与ART算法结合起来，用上一级网格单元的波慢给下一级网格单元的波慢赋值，重新计算每条射线穿过各个网格单元的射线长度，并利用ART算法计算和修改其波慢值，然后将该级网格继续划分，直到所有的网格单元均不能再划分为比成像单元更小的网格为止。塔式ART算法有效提高了计算的精度和图像重建质量，并有效反演出混凝土内部结构的强度分布以及缺陷的位置。

Description

混凝土超声层析成像算法

技术领域

本发明属于混凝土超声无损检测领域，特别涉及混凝土超声层析成像算法。

技术背景

层析成像(Computerized Tomography)是在不损伤研究“对象”内部结构的条件下，利用某种射线源，根据从“对象”外部用检测设备所获得的投影数据，依照一定的物理和数学关系，利用计算机反演“对象”内部未知的某种物理量的分布，生成二维、三维图像，重现“对象”内部特征。层析成像技术应用于混凝土无损检测的主要目标是在不损伤混凝土内部结构的条件下，确定建筑物内部的精细结构和局部不均匀性。

目前，比较常用的射线层析成像反演方法主要有反投影重建算法、变换重建类算法、离散图像代数重建类算法。其中，反投影算法是迄今为止速度最快的一种算法，但是计算精度不高，容易产生“伪迹”；以Fourier变换重建算法和滤波反(逆)投影算法为代表的变换重建方法，抗噪音干扰的能力差，并且如果投影数据不是沿直线的简单积分，那么可能就得不到解析反演公式的闭合形式，在这样的情况下，变换法就变得无效，因此不适用于混凝土层析成像；离散图像代数重建类算法中较为常用的有ART、SIRT、约束最小二乘类算法(含最大熵算法和最优化图像重建算法)等，适用于投影数据不完全，分布不均匀，射线路径为曲线的场合，并且便于计算机实现，因此得到广泛应用。其中，SIRT算法只有在测量数据特别不准确时才显示出它在重建质量上的优越性，而其他优点并不明显。而ART算法相对于约束最小二乘类算法，步骤简单，易于编程实现，因此，目前在混凝土超声层析成像中最常用的算法为ART算法。ART算法经过不断的改进，尽管仿真结果和试验结果是有效的，但图像重建的精度和快速性仍然没有很好的解决。

发明内容

本发明的目的在于克服上述现有技术不足，提出混凝土超声层析成像算法，本方法可以提高计算精度与计算速度，通过二维反演成像的结果，有效反映混凝土内部结构，从而确定缺陷的特征、大小和位置。

本发明的技术方案是这样实现的：混凝土超声层析成像算法，采用塔式ART算法按以下步骤进行：

步骤1：把所有的成像单元看作一个网格单元，作为第一级网格，根据经验赋予初始波慢

其中， 介于3500m/s与5000m/s之间；

步骤2：如果在当前级网格中存在某一网格单元G，其大小超过了成像单元大小，则按如下方法将该网格单元进行细分，并给新单元的慢度赋值：把一个成像单元看作一个像素，网格的长度和宽度都以像素为单位，并且都为整数，“[]”为向零取整的符号，网格单元G的大小表示为M×N个像素，如果网格单元G的长度和宽度均超过成像单元的长度与宽度时，该网格单元被近似平分成四个单元，按从左到右、从下到上的顺序其大小分别为[M/2]×[N/2]、(M-[M/2])×[N/2]、[M/2]×(N-[N/2])和(M-[M/2])×(N-[N/2])，并用该网格单元的波慢给这四个单元的波慢赋值；当网格单元的长度超过成像单元的长度，而其宽度与成像单元宽度相同时，该网格单元只在长度方向上近似平分成两个单元，按从左到右的顺序，其大小分别为[M/2]×1和(M-[M/2])×1，并用该网格单元的波慢给这两个单元的波慢赋值；当网格单元的宽度超过成像单元的宽度，而其长度与成像单元的长度相同时，该网格单元只在宽度方向上近似平分成两个单元，按从下到上的顺序，其大小分别为1×[N/2]和1×(N-[N/2])，并用该网格单元的波慢给这两个单元的波慢赋值，如果当前级的网格单元均不能再划分为比成像单元更小的网格时，则转到步骤5；

步骤3：对由步骤2得到的新一级的网格单元，按照直射线路径重新计算每条射线穿过各个网格单元的射线长度，即先确定各个网格单元所包含的成像单元，然后将各条射线穿过这些成像单元的长度相加，得到各条射线穿过各个网格单元的长度，于是，得到新的投影矩阵A；

步骤4：记第q轮迭代时第i条射线对第j个网格单元的波慢估算值为应用式(1)，逐条射线i(i＝1，2，L，n)逐轮对波慢作如下修改，其中0＜μ≤1，a_ij为步骤3求得的投影矩阵A的元素，每轮迭代完成后，判断是否满足收敛条件：‖f^q-f^q-1‖_∞＜ε，其中，

为第q轮迭代得到的慢度向量，ε为设定的误差界，是一个正数，若满足收敛条件，则停止本次迭代，转到步骤2，若不满足收敛条件，则继续按式(1)迭代。

${\hat{f}}_j^{q,i+1}={\hat{f}}_j^{q,i}+\mathrm{\μ}\frac{a_{\mathrm{ij}}}{\underset{j=1}{\overset{m}{\mathrm{\Σ}}}{a}_{\mathrm{ij}}^2}({\mathrm{\τ}}_i-{\widehat{\mathrm{\τ}}}_i^q)---\left(1\right)$

步骤5：停止计算，并输出每一成像单元的波速。

本发明可直接应用于混凝土无损检测现场，对被检测物体进行实时二维反演成像，从而准确反映混凝土内部结构。

本发明采用反演算法为塔式ART算法。将网格逐步细分，把上一级网格单元的波慢赋给下一级网格单元，重新计算每条射线穿过各个网格单元的射线长度，并利用ART算法计算和修改其波慢值，然后将该级网格继续划分，直到所有的网格单元均不能再划分为比成像单元更小的网格为止。通过不断细化网格，从而达到重建图像的目的。该算法运算速度快，成像效果比传统ART算法的成像效果好，缺陷位置更加准确和突出。

附图说明

图1是本发明塔式ART算法流程图；

图2是本发明单测检测方式示意图；

图3是本发明计算机模拟实验模型截面图；

图4是本发明计算机模拟实验的波速三维显示图，

其中，图(a)是模型波速三维显示图，图(b)是用传统的ART算法迭代100次得到的波速三维显示图，图(c)是用塔式ART算法计算得到的波速三维显示图；

图5是混凝土试件用传统ART算法迭代100次得到的断面波速分布图；

图6是混凝土试件用塔式ART算法计算得到的断面波速分布图。

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

具体实施方式

传统的ART算法中，网格划分一旦确定，投影矩阵A就惟一确定了，在迭代的过程中，每条射线穿过各个网格的长度固定不变；塔式ART算法将网格的动态划分与ART算法结合起来。每划分一次网格，重新计算每条射线穿过各个网格单元的长度，然后用ART算法计算并修正各个网格的波慢，直到网格不能再细分为止。

参照图1所示，其具体步骤如下：

1)把所有的成像单元看作一个网格单元，作为第一级网格，根据经验赋予初始波慢

其中，

介于3500m/s与5000m/s之间。

2)如果当前级的任一网格单元大小超过了成像单元大小，则将该级网格再细分，形成新一级的网格单元，并用上一级网格单元的波慢给新一级网格单元的波慢赋值；如果当前级的网格单元均不能再划分为比成像单元更小的网格时，则转到第5)步。

假设一个成像单元代表一个像素，网格的长度和宽度都以像素为单位，并且都为整数，符号“[]”代表向零取整。对大小为M×N个像素的网格单元细分的方法是：当M＞1时，则把该网格单元在长度方向上划分成2个单元，其长度分别为[M/2]和(M-[M/2])；当N＞1时，则把该网格单元在宽度方向上划分成2个单元，其宽度分别是[N/2]和(N-[N/2])。

当网格单元的长度和宽度均超过成像单元的长度与宽度时，该网格单元被近似平分成4个单元，按从左到右、从下到上的顺序其大小分别为[M/2]×[N/2]、(M-[M/2])×[N/2]、[M/2]×(N-[N/2])和(M-[M/2])×(N-[N/2])，并用该网格单元的波慢给这4个单元的波慢赋值；当网格单元的长度超过成像单元的长度，而其宽度与成像单元宽度相同时，该网格单元只在长度方向上近似平分成2个单元，按从左到右的顺序，其大小分别为[M/2]×1和(M-[M/2])×1，并用该网格单元的波慢给这2个单元的波慢赋值；当网格单元的宽度超过成像单元的宽度，而其长度与成像单元的长度相同时，该网格单元只在宽度方向上近似平分成2个单元，按从下到上的顺序，其大小分别为1×[N/2]和1×(N-[N/2])，并用该网格单元的波慢给这2个单元的波慢赋值。这样，对当前级的所有网格单元划分完毕后，就形成了下一级的网格单元，并且下一级的每个网格单元具有和其上一级父单元相同的波慢值。

3)对2)中得到的新一级的网格单元，按照直射线路径重新计算每条射线穿过各个网格单元的射线长度，即先确定各个网格单元所包含的成像单元，然后将各条射线穿过这些成像单元的长度相加，得到各条射线穿过各个网格单元的长度，于是，得到新的投影矩阵A。

4)记第q轮迭代时第i条射线对第j个网格单元的波慢估算值为 应用式(1)，逐条射线i(i＝1，2，L，n)逐轮对波慢作如下修改，其中0＜μ≤1。每轮迭代完成后，判断是否满足收敛条件：‖f^q-f^q-4‖∞＜ε，其中，

为第q轮迭代得到的慢度向量，ε为设定的误差界，是一个正数。若满足收敛条件，则停止本次迭代，转到第2)步，若不满足收敛条件，则继续按式(1)迭代。

${\hat{f}}_j^{q,i+1}={\hat{f}}_j^{q,i}+\mathrm{\μ}\frac{a_{\mathrm{ij}}}{\underset{j=1}{\overset{m}{\mathrm{\Σ}}}{a}_{\mathrm{ij}}^2}({\mathrm{\τ}}_i-{\widehat{\mathrm{\τ}}}_i^q)---\left(1\right)$

5)停止计算，并输出每一成像单元的波速。

本发明具有以下优点：塔式ART算法具有较快的迭代收敛速度及较高的计算精度，提高了图像重建质量，有效反演出混凝土内部结构的强度分布以及缺陷位置。

参照图2所示，在测区长度方向的两侧分别布置发射换能器(R₁，R₂，L R_m)以及接收换能器(T₁，T₂，L T_m)，每个换能器放在网格单元边界的中点处，其中，T_i表示第i个发射换能器，R_j表示第j个接收换能器。

本发明的实施例一：

参照图3所示，测区为100cm×60cm的长方形区域，划分为10×6的网格(成像单元)，其大小均为10cm×10cm，将网格按照从下到上、从左到右的顺序编号，其中22、23、28、29、34、35号为缺陷单元格，缺陷波速为4050m/s，正常波速为4500m/s，在长度方向的两侧布置发射以及接收换能器。

参照图4a-c所示，用ART算法迭代100次和塔式ART算法反演，计算结果对比如表1所示，图中，x轴代表测区长度，单位是米，y轴代表测区宽度，单位是米，z轴代表波速，单位是米/秒。可见，塔式ART算法在一定程度上削弱了异常体对周边单元波速的影响，与传统ART计算结果相比，缺陷单元周边地区的波速有所降低，缺陷单元波速更加接近真实值，计算精度有所提高。

                            表1塔式ART与传统ART计算结果对比

算法	正常区域波速			缺陷区域波速			波速离差	波速最坏距离
									最大值	最小值	平均值	最大值	最小值	平均值
	ART(100)	4622.46	4435.67	4489.47	4144.54	4124.5									4135.77	69.56	122.46
塔式ART							4616.85	4395.84	4492.47	4131.36	4091.62	4111.33	58.58	116.85

本发明的实施例二：

用混凝土试件做模拟试验，试件断面测区长44cm、深50cm，划分为11×10个网格，其中，中心部分为缺陷(混凝土夹泥)，大小为15×15cm，分别在长度方向的两侧布置发射以及接收换能器，获得11×11个声时，分别利用ART算法和塔式ART算法反演，得到的断面波速分布如附图5和附图6所示。

参照图5所示，传统ART算法的断面成像有所偏差，试件中心出现了两个缺陷区域，与混凝土试件中心为一均匀缺陷不尽相符。

参照图6所示，塔式ART算法成像质量有所提高，缺陷位置更加突出，但不能准确反映缺陷形状。

Claims (1)

1、混凝土超声层析成像算法，其特征在于，塔式ART算法按以下步骤进行：

步骤一：把所有的成像单元看作一个网格单元，作为第一级网格，根据经验赋予初始波慢 其中，

介于3500m/s与5000m/s之间；

步骤二：如果在当前级网格中存在某一网格单元G，其大小超过了成像单元大小，则按如下方法将该网格单元进行细分，并给新单元的慢度赋值：

把一个成像单元看作一个像素，网格的长度和宽度都以像素为单位，并且都为整数，式中的“[ ]”为向零取整的符号，网格单元G的大小表示为M×N个像素，如果网格单元G的长度和宽度均超过成像单元的长度与宽度时，该网格单元被近似平分成四个单元，按从左到右、从下到上的顺序其大小分别为[M/2]×[N/2]、(M-[M/2])×[N/2]、[M/2]×(N-[N/2])和(M-[M/2])×(N-[N/2])，并用该网格单元的波慢给这四个单元的波慢赋值；当网格单元的长度超过成像单元的长度，而其宽度与成像单元宽度相同时，该网格单元只在长度方向上近似平分成两个单元，按从左到右的顺序，其大小分别为[M/2]×1和(M-[M/2])×1，并用该网格单元的波慢给这两个单元的波慢赋值；当网格单元的宽度超过成像单元的宽度，而其长度与成像单元的长度相同时，该网格单元只在宽度方向上近似平分成两个单元，按从下到上的顺序，其大小分别为1×[N/2]和1×(N-[N/2])，并用该网格单元的波慢给这两个单元的波慢赋值，如果当前级的网格单元均不能再划分为比成像单元更小的网格时，则转到步骤五；

步骤三：对由步骤二得到的新一级的网格单元，按照直射线路径重新计算每条射线穿过各个网格单元的射线长度，即先确定各个网格单元所包含的成像单元，然后将各条射线穿过这些成像单元的长度相加，得到各条射线穿过各个网格单元的长度，于是，得到新的投影矩阵A；

步骤四：记第q轮迭代时第i条射线对第j个网格单元的波慢估算值为

应用式(1)，逐条射线i(i＝1，2，L，n)逐轮对波慢作如下修改，其中0＜μ≤1，a_ij为步骤3求得的投影矩阵A的元素，每轮迭代完成后，判断是否满足收敛条件：‖f^q-f^q-1‖∞＜ε，其中，

为第q轮迭代得到的慢度向量，ε为设定的误差界，是一个正数，若满足收敛条件，则停止本次迭代，转到步骤二，若不满足收敛条件，则继续按式(1)迭代，

${\hat{f}}_j^{q,i+1}={\hat{f}}_j^{q,i}+\mathrm{\μ}\frac{a_{\mathrm{ij}}}{\underset{j=1}{\overset{m}{\mathrm{\Σ}}}{a}_{\mathrm{ij}}^2}({\mathrm{\τ}}_i-{\widehat{\mathrm{\τ}}}_i^q)---\left(1\right)$

步骤五：停止计算，并输出每一成像单元的波速。

Images