CN104931584A - 基于压缩采样理论的超声ct检测方法 - Google Patents

Abstract

本发明涉及一种基于压缩采样理论的超声CT检测方法，该方法采用压缩采样方法进行优化，改进了传统的超声CT检测方法，可减少数据的采集量。其实施方法为：首先用计算机生成一个随机的测量矩阵，根据测量矩阵确定测量路径，并将发射探头和接收探头分别放置在试件两侧，一端发射一端接收，得到各路径的传播时间，将各个路径的传播时间按照对应的顺序带入到l1优化算法中，最终对混凝土试件内部损伤情况进行重构。本发明所设计的基于压缩采样的超声CT方法，只需要少量的测量就可以高精度地重构出试件内部的结构，找到损伤位置并且确定损伤的大小，能够大大地加快损伤定位的速度，提高了监测效率，减少了工作量，具有良好的应用前景。

Description

基于压缩采样理论的超声CT检测方法

技术领域

本发明专利涉及一种基于压缩采样理论的超声CT方法。

背景技术

压缩采样是由E.J.Candès、J.Romberg、T.Tao和D.L.Donoho等科学家于2004年提出的，其主要应用于图像数据的获取和压缩，压缩采样就是在数据采样时直接采集压缩格式的数据，而不是先采集全部的数据然后再进行压缩处理。这样可以减少数据获取的成本，有利于数据的传输与保存。压缩采样与Nyquist采样定理相比较，采样定理是在模态/数字信号的转换过程中，当采样频率大于信号最高频率的2倍时候，采样之后的数字信号才能完整地保留原始信号的信息。如果信号频率较高的话，那么采样得到的数据量也是巨大的。而压缩采样和采样频率无关，只和信号的稀疏性有关，如果信号是稀疏的或者经过某种变换之后是稀疏的，那么只需要随机采样即可以保留原始信号中的信息，进而通过算法对原始信号进行恢复。

在土木结构健康监测中，传统的超声CT法尽管具有人为误差小，监测准确等优点，但是这种传统的方法会耗费极大的人力物力，效率低下，实际应用性差。因此，在建筑行业飞速发展的今天，更多重要的枢纽性建筑如主梁，柱，大坝等迫切需要一种高效的监测方法来保证其在运营期间的安全。

发明内容

基于以上不足之处，本发明公开一种基于压缩采样理论的超声CT方法，应用于监测高铁无砟轨道，具体步骤如下：

第一步：确定测量精度，将要监测的试件在监测平面划分成数个单元；

第二步：用计算机生成一个随机的测量矩阵Φ，用于确定测量路线，并用Matlab计算出由各测量路线穿过试件所经过的各单元长度构成字典矩阵A；

第三步：根据上步得到的测量矩阵，从传统的超声CT方法的全测量路径中随机地挑选出本次测量路径；

第四步：按照第三步得出的测量路径在试件两侧放置发射器和接收器，得到测量路径的传播时间

第五步：重复第四步，得到所有选择路径的传播时间，进而构成此次检测的测量向量T′_m；

第六步：建立试件重构模型如下式：

ΔT′_m＝T′_m-T₀＝Φ(A·S-A·S₀)＝Φ·A·ΔS   (1)

其中，ΔT′_m为实际检测的传播时间向量T′_m与相应路径构成的无损伤状态下的标准传播时间向量T的差向量，被检测试件的慢度向量s，其与无损伤状态下的标准慢度向量的差为慢度差向量ΔS，由于损伤在试件中是稀疏分布的，因此该慢度向量ΔS也是稀疏的向量；

第七步：采用基于压缩采样的凸优化方法进行求解，表达式如下

$\mathrm{\Δ}\hat{S}=\arg \min ({\left|\right|\mathrm{\Θ}\·\mathrm{\ΔS}-\mathrm{\ΔT}\left|\right|}_2+\mathrm{\λ}{\left|\right|\mathrm{\ΔS}\left|\right|}_1)---\left(2\right)$

进而有

$\hat{S}=S_0+\mathrm{\Δ}\hat{S}---\left(3\right)$

根据重构出的慢度向量还原试件内部损伤情况。

本发明还具有如下技术特征：

进行网格划分时，将监测区划分为p×q＝w个规则网格，w是网格的总数，其中p、q为划分的矩形单元网格数，监测区共有N条射线通过，L_i为第i条射线，T_i是第i条射线的旅行时间，定义速度的倒数为慢度s＝1/V，由Radon公式：

$T_i={\∫}_{L_i}\frac{1}{V_j(x,y)}\mathrm{dL}={\∫}_{L_i}{s}_j(x,y)\mathrm{dL}---\left(4\right)$

其中：

t_ij——是第i条射线在第j个单元中传播的时间

V_j(x，y)——第j个单元中波的传播速度

s_j(x，y)——第j个单元中波的传播的慢度

L_i——第i条射线的长度

设网格足够细，成像单元足够小，每一个单元中的速度V_j(x，y)和慢度s_j(x，y)均视为为一个常数，将上式的积分写成如下累计求和的形式：

$T_i=\underset{j=1}{\overset{w}{\mathrm{\Σ}}}{t}_{\mathrm{ij}}=\underset{j=1}{\overset{w}{\mathrm{\Σ}}}\frac{a_{\mathrm{ij}}}{V_j}=\underset{j=1}{\overset{w}{\mathrm{\Σ}}}{a}_{\mathrm{ij}}{s}_j---\left(5\right)$

探测中，总共发出的N条超声波射线，每一条射线到达对面的行进时间构成了一个时间向量T＝(T₁ T₂…T_n)^T，形成成如下方程组的形式：

$\left\{\begin{array}{c}{T}_1=\underset{j=1}{\overset{w}{\mathrm{\Σ}}}\frac{a_{1j}}{V_j}=\frac{a_{11}}{V_1}+\frac{a_{12}}{V_2}+...+\frac{a_{1w}}{V_w}\\ T_2=\underset{j=1}{\overset{w}{\mathrm{\Σ}}}\frac{a_{2j}}{V_j}=\frac{a_{21}}{V_1}+\frac{a_{22}}{V_2}+...+\frac{a_{2w}}{V_w}\\ .\\ .\\ .\\ T_n=\underset{j=1}{\overset{w}{\mathrm{\Σ}}}\frac{a_{\mathrm{nj}}}{V_j}=\frac{a_{n1}}{V_1}+\frac{a_{n2}}{V_2}+...+\frac{a_{\mathrm{nw}}}{V_w}\end{array}\right.---\left(6\right)$

改写成慢度的形式：

$\left\{\begin{array}{c}{T}_1=\underset{j=1}{\overset{w}{\mathrm{\Σ}}}{a}_{1j}{s}_j=a_{11}{s}_1+a_{12}{s}_1+...+a_{1w}{s}_w\\ T_2=\underset{j=1}{\overset{w}{\mathrm{\Σ}}}{a}_{2j}{s}_j=a_{21}{s}_1+a_{22}{s}_1+...+a_{2w}{s}_w\\ .\\ .\\ .\\ T_n=\underset{j=1}{\overset{w}{\mathrm{\Σ}}}{a}_{\mathrm{nj}}{s}_j=a_{n1}{s}_1+a_{n2}{s}_1+...+a_{\mathrm{nw}}{s}_w\end{array}\right.---\left(7\right)$

改写成矩阵的形式：

T＝AS   (8)

其中：

T——是所有N次检测超声传播的时间向量；

A——是所有N次检测在w个单元中的长度矩阵；

a_ij——超声波在第i个路径上经过第j个单元的长度；

S——所有w个单元的慢度向量。

随机测量路径的选择由测量矩阵Φ来决定，随机测量矩阵Φ中的所有元素都为0或1，并且每一行只有一个元素为1其余元素均为0，如果第r行第j列的元素为1，意味着选取第j条路径作为第r次随机观测的观测路径，根据测量矩阵Φ从传统的N条全测量路线中随机选取取m条测量路线，其中m满足压缩采样的关系m＞μ·K·log(N/K)，μ因问题不同取值不同，μ≈4，这样试件的内部情况就能够被精确重构。在测量中，每个路径读三个传播时间，取平均值作为该路径信号的传播时间。

本发明的特点和优点：

本发明克服了传统超声CT检测方法具有的弊端即需要对试件进行大量冗余检测的缺点。该方法只需要少量的测量就能以较高精度重构出试件内部的结构，找到损伤位置和大小，大大地加快了损伤定位的速度，极大地提高监测效率。

附图说明

图1为灌装柱的超声CT成像结果和示意图

图2为有孔洞损伤结构的超声波的绕行路径

图3为超声波传播路径经过裂缝缺陷的一种典型情况

图4为由计算机生成的一次随机测量矩阵Φ

图5为压缩采样表达式图解

图6为基于压缩采样的超声CT法流程图

图7为试件设计主视图

图8为试件设计侧视图

图9为激励信号波形

图10为超声波经过不同路径传播的波形

图11为由测量矩阵得到的测量路线

图12为压缩采样方法的重构结果

图13为试件情况及网格划分

图14为测量矩阵和随机选取的测量路径

图15为情况1无损情况下的试件的识别情况

图16为情况2试件损伤情况

图17为情况2试件的实际慢度矩阵与无损慢度矩阵的差

图18为情况2有损与无损试件慢度矩阵差的识别结果

图19为情况3试件损伤情况

图20为情况3有损与无损试件慢度矩阵差的识别结果

具体实施方式

实施例1

当超声波或声波通过物体传播时，物体内部的单元会与波动信号相互作用，对声波的速度产生影响。如果物体内部某处比较致密且坚硬则波动信号传播的速度就会很快，反之，则会减慢，某条测线的波速实际上是构成该条测线各部份波速的综合值，可以计算出声波在物体中传播的速度，来确定物体内部损伤的情况。将物体层(断)面划分成一定数量的网格(亦称成像单元)，一侧激发，另一侧所有点接收，各成像单元被测线多次穿过。采用迭代方法反演各成像单元的波速值，重建物体内部波速图像，这就是声波层析成象的原理及实施技术。图1a)所示是在一段灌注桩中进行CT扫描测量。图1b)是由扫描获得的断面波速等值线色谱图，图中的低波速区正是被包裹在混凝土中的泥砂团。

将监测区划分为p×q＝w个规则网格(其中p、q为划分的矩形单元网格数，w是网格的总数)，假设测区共有N条射线(超声波行进路径)通过，L_i为第i条射线，T_i是第i条射线的旅行时间，定义速度的倒数为慢度s＝1/V，由Radon公式：

$T_i={\∫}_{L_i}\frac{1}{V_j(x,y)}\mathrm{dL}={\∫}_{L_i}{s}_j(x,y)\mathrm{dL}---\left(1\right)$

其中：

t_ij——是第i条射线在第j个单元中传播的时间

V_j(x，y)——第j个单元中波的传播速度

s_j(x，y)——第j个单元中波的传播的慢度

L_i——第i条射线的长度

假设网格足够细，成像单元足够小，每一个单元中的速度V_j(x，y)和慢度s_j(x，y)均视为为一个常数。则可将上式的积分写成如下累计求和的形式：

$T_i=\underset{j=1}{\overset{w}{\mathrm{\Σ}}}{t}_{\mathrm{ij}}=\underset{j=1}{\overset{w}{\mathrm{\Σ}}}\frac{a_{\mathrm{ij}}}{V_j}=\underset{j=1}{\overset{w}{\mathrm{\Σ}}}{a}_{\mathrm{ij}}{s}_j---\left(2\right)$

探测中，总共发出的N条超声波射线，每一条射线到达对面的行进时间构成了一个时间向量T＝(T₁ T₂…T_n)^T，可以写成如下方程组的形式。

$\left\{\begin{array}{c}{T}_1=\underset{j=1}{\overset{w}{\mathrm{\Σ}}}\frac{a_{1j}}{V_j}=\frac{a_{11}}{V_1}+\frac{a_{12}}{V_2}+...+\frac{a_{1w}}{V_w}\\ T_2=\underset{j=1}{\overset{w}{\mathrm{\Σ}}}\frac{a_{2j}}{V_j}=\frac{a_{21}}{V_1}+\frac{a_{22}}{V_2}+...+\frac{a_{2w}}{V_w}\\ .\\ .\\ .\\ T_n=\underset{j=1}{\overset{w}{\mathrm{\Σ}}}\frac{a_{\mathrm{nj}}}{V_j}=\frac{a_{n1}}{V_1}+\frac{a_{n2}}{V_2}+...+\frac{a_{\mathrm{nw}}}{V_w}\end{array}\right.---\left(3\right)$

改写成慢度的形式：

$\left\{\begin{array}{c}{T}_1=\underset{j=1}{\overset{w}{\mathrm{\Σ}}}{a}_{1j}{s}_j=a_{11}{s}_1+a_{12}{s}_1+...+a_{1w}{s}_w\\ T_2=\underset{j=1}{\overset{w}{\mathrm{\Σ}}}{a}_{2j}{s}_j=a_{21}{s}_1+a_{22}{s}_1+...+a_{2w}{s}_w\\ .\\ .\\ .\\ T_n=\underset{j=1}{\overset{w}{\mathrm{\Σ}}}{a}_{\mathrm{nj}}{s}_j=a_{n1}{s}_1+a_{n2}{s}_1+...+a_{\mathrm{nw}}{s}_w\end{array}\right.---\left(4\right)$

改写成矩阵的形式：

T＝AS   (5)

其中：

T——是所有N次检测超声传播的时间向量

A——是所有N次检测在w个单元中的长度矩阵(元素各不相同，变化很大)

a_ij——超声波在第i个路径上经过第j个单元的长度

S——所有w个单元的慢度向量(速度的倒数)

目前工程CT大多采用迭代方法，在层析成像问题中，用得较多的迭代方法有代数重建技术(ART)、最优化准则、联合迭代技术(SIRT)和阻尼最小二乘法(LSQR)。或者非迭代方法——反投影技术(BPT)生成迭代算法所需的波慢度初始值等来计算，以此来重构试件内部的损伤情况。

超声波在混凝土中的传播速度为3,800-4,500m/s，而在钢筋和空气中则分别以5,200m/s和340m/s的速度传播。考虑到超声波在空气中能量大幅度衰减，传播过程中遇到的孔洞裂缝这一类缺陷应该被视为障碍。然而，由于钢筋与混凝土都是固体，超声波在二者之中传播速度比较相似。因此，超声波从混凝土中射入钢筋，再由钢筋射入混凝土。鉴于超声波在钢筋中的传播速度快于混凝土，在钢筋中的传播时间较于混凝土要小一些，我们可以用这种现象来鉴别混凝土中的钢筋。基于惠更斯-菲涅尔原理，当超声波遇到阻碍或裂缝时会发生衍射，在遇到阻碍时会停止传播，障碍的边缘会变成一个新的球形波源，这些新的球面波决定了最后接收时刻的波形。换句话说，波为了穿过缺陷或障碍，只能在它们的边缘传播。

如图2a)所示，在传感器的发射端到接收端的连线(虚线)上有一个障碍，传感器发射端发出的超声波经过障碍的边缘传播到传感器接收端，如图2a)中的实线所示。虚线是用于确定标准传播时间T_{0_i}的标准路径，而实线的长短则决定真实传播路径的传播时间。

如图2b)所示，当传感器发射端和接收端之间连线上存在两个障碍时，超声波从第一个障碍边缘经过并发生了第一次衍射，接着由于直接传播路径穿过障碍发生了第二次衍射。最终，超声波的实际传播路径如实线所示。因此，当传播路径中存在障碍或裂缝时，传播时间随传播路径长度增加而变长。总的来说，传播时间的变化量ΔT_i和直接传播路径中障碍和裂缝的数量存在着叠加关系。

此外，考虑到损伤的多种可能种类，图3给出了损伤种类是裂缝时，损伤对测量路径的影响情况，两个典型路径分别为路径1(接收探头2和发射探头2所连平行于裂缝的虚线)和路径2(接收探头1和发射探头1所连垂直于裂缝的虚线)。值得注意的是，孔洞可以被看做混凝土结构中的各向同性损伤，即不同测量路径的选择对传播时间变化量影响不大。然而，裂缝却是各向异性的缺陷，裂缝对路径1和路径2影响程度如图3所示，可见缺陷的存在对传播路径的影响程度有很大不同。据此，在对裂缝进行定位之后，可以通过观察时间变化量最小ΔT_min的路径来确定裂缝方向。另外，传播时间变化量最大ΔT_max的路径垂直于裂缝方向。因此，我们可以从裂缝中分辨出孔洞进而确定裂缝的方向。

压缩采样是近年来最新发展的数据压缩技术，由美国加州理工学院的教授等于2005年提出，目前已经成为热门的研究方向。该方法主要应用于图像数据的获取和压缩，所谓的压缩采样就是在数据采样时直接采集压缩格式的数据，而不是先采集全部的数据然后再进行压缩处理。这样可以减少数据获取的成本，有利于数据的传输与保存。压缩采样与Nyquist采样定理相比较，采样定理是在模态/数字信号的转换过程中，当采样频率大于信号最高频率的2倍时候，采样之后的数字信号才能完整地保留原始信号的信息。如果信号频率较高的话，那么采样得到的数据量也是巨大的。而压缩采样和采样频率无关，只和信号的稀疏性有关，如果信号是稀疏的或者经过某种变换之后是稀疏的，那么只需要随机采样即可以保留原始信号中的信息。

信号的稀疏性是压缩采样的一个重要特征。如果一个一维的向量X其大小为N×1属于实数域其元素中只有K个非零的量，则称向量X为K-稀疏的向量；如果信号本身不稀疏，但是通过某种变换或者在某个基矩阵亦或字典矩阵的作用下是稀疏的，如下式所示，则信号依然是可压缩的。信号越稀疏，则其可压缩性就越好。

α＝Ψ′x，x＝Ψα  (6)

其中，

α——信号的分解系数，其为稀疏向量

Ψ——基矩阵，可以为某种变换的正交基矩阵，也可以是字典矩阵

Ψ′——基矩阵的逆矩阵，如果基矩阵为正交基矩阵则其逆矩阵Ψ′＝Ψ^T假设信号x通过线性映射的方式，得到长度为m(m＜N)的向量y。

y＝Φx

y＝ΦΨα＝Θα   (7)

其中，

Φ——m×N的测量矩阵或者感知矩阵，需具有一定的随机性。

Θ——m×N变换矩阵

一般来说，因为m＜N，所以上式是病态问题，无法直接从y重构出X。但是如果X在矩阵Ψ下的系数向量α是稀疏或者近似稀疏的则上述问题可以通过求解凸优化的问题重构出α，进而恢复信号X。即在满足一定条件下，可以通过求解如下带有约束条件的线性规划问题，从y中精确恢复α。

$\begin{array}{ccc}\widehat{\mathrm{\α}}=\min {\left|\right|\mathrm{\α}\left|\right|}_1& \mathrm{subject\; to}& y=\mathrm{\Θ\α}---\left(8\right)\end{array}-$

其中是通过优化算法求得的α的唯一精确解。

经研究发现，变换矩阵Θ必须要满足约束等距特性RIP(RestrictedIsometry Property)的条件，稀疏向量α才能通过优化算法得以精确重构。

矩阵Θ的约束等距特性是指对于整数s∈{1，2，…，N}，矩阵Θ存在一个等距常数δ_s，这里δ_s是满足下式的一个最小的数，且对于所有稀疏度为K的向量v_s都成立。

$1-{\mathrm{\δ}}_s\≤\frac{{\left|\right|{\mathrm{\Θv}}_s\left|\right|}_2}{{\left|\right|v_s\left|\right|}_2}\≤1+{\mathrm{\δ}}_s,{\mathrm{\δ}}_s>0---\left(9\right)$

如果矩阵Θ满足RIP条件且α是稀疏的，则是通过优化算法求得的α的唯一精确解，x也可以被精确重构。要寻求满足条件的测量矩阵Φ，首先在基矩阵Ψ不变的情况下，需要构造矩阵，使得Θ＝ΦΨ满足RIP条件。此外，需要寻找具有任意选取列向量构成的子集都具有近似正交性的感知矩阵，且这些子集越大越好，也就是测量矩阵Φ与基矩阵Ψ的相干性越小越好。然而，单纯采用不相干性和RIP条件设计观测矩阵是一个病态问题，在实际上基本不可行。幸运的是，研究表明随机观测矩阵能以很大的概率同时满足不相干性和RIP条件，只要满足m＞μ·K·log(N/K)，其中μ随问题不同取值不同，一般μ≈4。故实际中一般采用计算机生成随机矩阵进行观测，常见的有高斯观测矩阵、伯努利矩阵等。

在实际应用中，由于环境因素等的影响，测量信号难免被噪声所污染，因此，在此种情况下，考虑噪声的影响，重构方程变为：

$\begin{array}{ccc}\widehat{\mathrm{\&alpha;}}=\arg \min {\left|\right|\mathrm{\&alpha;}\left|\right|}_1& \mathrm{subject\; to}& {\left|\right|\mathrm{\&Theta;\&alpha;}-y\left|\right|}_2<\mathrm{\&epsiv;}---\left(10\right)\end{array}$

本发明主要是针对传统超声CT法检测混凝土试件需要对被测试件进行大量的，冗余的检测这一问题，并结合压缩采样理论，提出一种应用于混凝土试件检测的基于压缩采样的超声CT方法。该方法只需要少量的测量就可以精度很高地重构出试件内部的结构，找到损伤位置和大小，大大地加快了损伤定位的速度，极大地提高了监测效率。

要监测混凝土试件内部情况，首先将要监测的混凝土试件划分成很多单元。令混凝土单元的慢度为标准慢度S₀，可知划分网格后的大部分的单元其慢度都是S₀，只有少量包含有缺陷的单元，其慢度不为S₀。如果单元的慢度大于S₀，则该单元有可能包含有孔洞或者裂缝；反之，该单元有可能有钢筋穿过。整个结构的所有单元的真实慢度矩阵S与标准慢度矩阵S₀做差如式(11)，那么差向量ΔS中多数元素为0，只有少数的元素不为0，符合压缩采样理论中的稀疏性的定义。

ΔS＝S-S₀   (11)

由于T＝A·S上式变为如下形式：

T＝A·S＝A(S₀+ΔS)＝A·S₀+A·ΔS＝T₀+ΔT   (12)

进而得到：

ΔT＝T-T₀＝A(S-S₀)＝A(ΔS)   (13)

式中的测量到的时间差向量ΔT和矩阵A分别对应于压缩采样理论的原始信号x和正交基矩阵Ψ。应该注意的是，传播距离矩阵A应该被视为字典矩阵，它并不满足正交性条件，然而根据前文所述理论，式(10)在A不满足正交性条件的情况下仍可以得出唯一解。式(13)可以写成下面形式

其中ΔT_i是超声波沿i路径的传播时间差，T_i和T_{0_i}分别是有损情况和无损情况的传播时间，a_ij是超声波在第i个路径上经过第j个单元的长度，Δs_j是第j个单元的慢度差。

在测量之前，预先用计算机生成一个用于确定测量路径的测量矩阵Φ，如图4所示，从N条全测量路径中通过测量矩阵Φ来选取m条测量路径(m＜N)，公示表达如下

ΔT′_m＝ΦΔT

$\mathrm{\&Phi;}={\left[\begin{array}{cccccc}0& & & 1& & \\ & & 1& & & \\ & & & & 1& \\ & 1& & & & \\ & & & & & 1\end{array}\right]}_{m\&times;N}---\left(15\right)$

其中，对角矩阵Φ中的元素非0即1，这些单元的分布由随机伯努利分布确定。根据有数值的单元位置相关的路径会被扫描，否则，路径将不会被扫描。例如，如果那么路径i会被扫描，如果那么路径j将不会被扫描。换句话说，测量矩阵中非零单元的存在与选出的路径是一致的。

将(13)式带入(15)式，时间向量可以表示如下：

ΔT′_m＝Φ·ΔT＝Φ·A·ΔS＝ΘΔS   (16)

其中字典矩阵A可由Matlab计算得出，以上过程可用图5表示。

对于测量矩阵Φ，应由计算机随机生成，首先使其与字典矩阵相乘而得到的变换矩阵Θ满足RIP条件，其次应该使测量矩阵与字典矩阵相关性尽可能小，但利用这两个条件设计观测矩阵是一个病态问题，在实际上基本不可行。幸运的是，研究表明随机观测矩阵能以很大的概率同时满足不相干性和RIP条件，只要满足m＞μ·K·log(N/K)，(μ≈4)。故实际中一般采用计算机生成随机矩阵进行观测，常见的有高斯观测矩阵、伯努利矩阵等。这样就可以从传统的超声CT法中的全测量路径中选取少量路径，进而很好地完成测量和重构。

ΔT′_m需要由有损情况测量路径传播时间向量T′_m和无损情况的传播时间向量T₀求得。其中T′_m可由少量m次测量得到，应按照测量矩阵规定的测量路径在试件两侧放置发射器和接收器，每个路径读三个传播时间，取平均值作为该路径信号的传播时间，并完成所有路径的测量得到向量T′_m。由于T′_m是一个长度为m的向量且长度小于慢度向量长度N，所以式(17)是病态问题，无法直接从T′_m重构出ΔS。但是向量ΔS满足稀疏性条件，且变换矩阵Θ满足RIP条件，则可以将测量得到的各路径的传播时间采用基于压缩采样的凸优化方法进行求解，先将测量得到的结果T_m减去测量的各路径的基本时间向量T₀＝d_i/V₀，得到测量值与无损情况的差信号时间向量ΔT′_m，再将差信号时间向量带入到压缩采样的重构算法(13)中进行求解，表达式如下

$\mathrm{\&Delta;}\hat{S}=\arg \min ({\left|\right|\mathrm{\&Theta;}\&CenterDot;\mathrm{\&Delta;S}-\mathrm{\&Delta;}{T}_m^{\&prime;}\left|\right|}_2+\mathrm{\&lambda;}{\left|\right|\mathrm{\&Delta;S}\left|\right|}_1),\mathrm{\&Theta;}=\mathrm{\&Phi;}\&CenterDot;A---\left(17\right)$

其中λ是l1最小化算法中的拉格朗日乘子。

求得的可以用来定位损伤位置，其中是通过优化算法求得的的唯一精确解，也是缺陷位置的标志。在中非的0单元的值与缺陷的大小和种类相关。

$\left|\mathrm{\&Delta;}\hat{S}\right|\&GreaterEqual;\mathrm{\&Delta;}{\stackrel{\&OverBar;}{S}}_0---\left(18\right)$

$\hat{S}=S_0+\mathrm{\&Delta;}\hat{S}---\left(19\right)$

是损伤定位的门槛值，式(18)意义即为大于这个临界值的时候可以将慢度差异识别为损伤。因此上述问题可以通过求解凸优化的问题重构出ΔS。即在满足一定条件下，可以通过求解如下带有约束条件的线性规划问题，从而精确恢复ΔS，这样试件的内部损伤就能够被精确重构。

值得注意的是，此算法需要远小于传统的超声CT检测方法的测量路径就可以对混凝土试件进行精确重构。

综合以上技术，基于压缩采样的超声CT法流程图如图6所示，其步骤可总结为：

第一步：确定测量精度，将要监测的试件在监测平面划分成数个单元。

第二步：用计算机生成一个随机的测量矩阵Φ，用于确定测量路线。并用Matlab计算出由各测量路线穿过试件所经过的各单元长度构成字典矩阵A。

第三步：根据上步得到的测量矩阵，从传统的超声CT方法的全测量路径中随机地挑选出本次测量路径。

第四步：按照第三步得出的测量路径在试件两侧放置发射器和接收器，一端发射一端接收，得到测量路径的传播时间

第五步：重复第四步，得到所有选择路径的路径的传播时间，进而构成此次检测的测量向量T′_m。

第六步：建立试件重构模型ΔT′_m＝T′_m-T₀＝Φ(A·S-A·S₀)＝Φ·A·ΔS

其中，ΔT′_m为实际检测的传播时间向量T′_m与相应路径构成的无损伤状态下的标准传播时间向量T的差向量。被检测试件的慢度向量s，其与无损伤状态下的标准慢度向量的差为慢度差向量ΔS。由于损伤在试件中是稀疏分布的，因此该慢度向量ΔS也是稀疏的向量。

第七步：采用基于压缩采样的凸优化方法进行求解，根据重构出的慢度向量还原试件内部损伤情况。

上述基于压缩采样理论的超声CT检测方法，可应用于混泥土试件如梁、柱、墙、桩基、大坝等，也可应用于金属构件的检测如铁、铝制构件。

实施例2

预设缺陷的混凝土试件设计图如图7、图8所示，试件采用C40混凝土浇筑。试块中分别用三个孔洞直径10mm，直径15mm，直径20mm的空心PVC管来模拟不同大小的孔洞，用三个等大的可抽离的木板(试块成型后抽离)分别模拟0°、45°和90°三个不同方向的裂缝，并且在不同深度埋置了2根直径为12mm的钢筋，3根直径为22mm的钢筋的带肋钢筋。试件的整体尺寸是1000mm(L)×300mm(T)×200mm(W)，将其横截面1000mm(L)×300mm(T)作为试验的探测区域。

本次试验采用的超声探测仪是北京宇通时代科技公司生产的NU62型非金属超声波检测仪，配备有两个收发两用的超声波探头，探头直径为4cm。声时精度为0.05μs，试验中的采样频率为20MHz，波形的主频率为50kHz，信号电压为250V，波形如图9所示。

本实验中探测精度取为2.5cm，值得一提的是，使用单通道同时探测所有路径的传播时间即是不需要的也是不可能实现的。幸运的是，重构的试件是静态物体，只需要在不同时间采用不同路径探测就可以解决这一问题。实验探测方法为：首先固定一端(比如底端)的驱动器探头，然后在对面的一端(上端)每隔5cm移动另外的一个作为传感器的超声探头，在仪器上读取两个点之间的时间延迟，读取三个数取平均值作为该组的传播时间；待本轮测量结束后，将底端的驱动器移动到下一个位置，重复上面的步骤，直到完成整个区域的探测。

试验中的纵向探测精度和横向探测精度分别为2.5cm和1cm，划分的单元总数为30×33＝900个。此外，在实验中考虑到试件的边缘不适宜放置传感器，并且会对传播时间起到一定的干扰作用，因此左右两端各留出了10cm的空白区域不予放置传感器探测。最终长度方向的探测距离是从10-90cm，总测量量为33×33＝1089个点，然而，基于压缩采样方法，只需采用200条测量路径去重构试件的内部情况已经足够，由测量矩阵得出的测量路径如图11所示，再根据测量路线，分别在试件两侧贴上发射器和接收器头。

从图10中可见，首波时间决定于传播路径长度和试件的内部条件，四个路径均为探头在时间两侧对正，且连线垂直于侧面。四条路径分别代表了超声波在时间中传播的四种最简单的情况。a路径为不经过任何损伤的路径；b路径为经过了一个圆孔状损伤的路径，可以看到超声波的绕行路线；c路径经过了一个竖直裂缝的路径，超声波沿着裂缝的长度方向绕行；d路径是超声波经过横向裂缝的路径，此时超声波取两条d-1和d-2路径绕走裂缝，从路径的长短可以大致发现d-1要比d-2路径短一些。从图中可以精确看出不同路径的首波时间。超声波沿无损的a路径的传播时间为70.72μs；经过有孔洞损伤的b路径传播时间为71.64μs；超声波经过竖向裂缝绕行c路径的传播时间为71.16μs；观察d路径，d-1路径传播时间为71.16μs和d-2路径传播时间为74.24μs。通过对比一系列路径传播时间，可见缺陷对首波时间的影响是很明显的，基于这种特点可以通过测量不同路径的首波时间来判定缺陷的类型。

完成所有需要测量的传播路径的传播时间的测量，得到测量路径的数据后，可以通过l1最小化方法得到有损情况试件慢度矩阵与无损矩阵的差值，且通过设置损伤识别门槛值，如式(18)所示，是损伤定位的门槛值，即大于这个临界值的时候可以将慢度差异识别为损伤，设置代表损伤直径小于1mm的损伤会被忽略掉，损伤识别的效果图如图12所示，可以看出钢筋对识别效果几乎没有影响，原因在于钢筋直径和混凝土中的缺陷相比可以忽略，且超声波在钢筋中的传播速度与混凝土无太大差异，最后我们定位了缺陷的位置和大小，它们分别为(200，100，9.8)，(300，200，15.2)，(400，100，20.5)和(650，100，19.4)还有三个明显的裂缝缺陷，这些都可以从图中清楚看出来。同样在识别结果中由于噪音，测量路线数目，选择测量路径数目等因素影响，都会造成损伤识别存在一些不可避免的偏差。

基于本专利方法，我们设计了三种不同的实验模式进行数值研究。情况1是完好无损的混凝土试件；情况2是有5个损伤的混凝土试件，试件的种类均为孔洞；情况3是有6个损伤的混凝土试件，损伤种类包括裂缝和孔洞。

我们采用一个基于二维混凝土结构平面的数值研究来证实本专利提出的算法。实验中采用一个长度100cm宽度30cm的混凝土厚板，采用的网格尺寸大小为25mm×10mm，共有w＝40×30＝1200个网格，如图13所示，Φ6150表示混凝土试件中的箍筋直径为6mm，箍筋间距为150mm。传统超声探测方法中，共有N＝40×40＝1600条探测路径。因此，必须指出的是，传统探测方法得到的数据有相当一部分是冗余的。

采用本专利提出的算法，从1600条测量路径中选取的200条测量路径，即m＝200，满足该算法要求的条件m＞μ·K·log(N/K)。因此，基于上述网格划分，测量矩阵Φ的维度为200行×1600列，那么变换矩阵Θ＝ΦA的维度为200行×1200列。变换矩阵满足不相干性和RIP条件。实验中，测量路径由随机测量矩阵确定并选取，测量矩阵如图14所示。

超声波以标准速度在混凝土中传播。如果传播路径中存在一个损伤，那么超声波会波及到附近的网格单元进而引起超声波传播时间的延长。

情况1：

在混凝土试件内部无损的情况下，超声波经过第i个测量路径的传播时间T_i与标准路径传播时间T_{0_i}相同，这意味着时间差向量ΔT是一个零向量。并且，由l1最小化算法计算得到的慢度差向量ΔS也是一个零向量。识别结果结果如图15所示，探测区域整片都是绿色的表示混凝土结构是健康的。

情况2：

我们在混凝土结构中设置了5个孔洞单元，其中有4个直径为2mm，其中一个孔洞3mm，且它们的位置分别为(100，150)，(200，230)，(350，100)，(750，50)和(550，250)。需要指出的是，网格尺寸大于孔洞尺寸，所以两种直径对应的慢度分别为1.0160和1.0240。实际结构的慢度矩阵与无损结构的慢度矩阵作差如图17所示。

如上所述，如图17所示超声波在遇到损伤的时候会绕过损伤继续传播。因此，在模拟中，我们根据惠更斯-菲涅尔原理重新计算了测量路径中有缺陷存在时的传播路径长度。例如：当超声波传播绕过大孔洞时，路径的长度为300.24mm，经过计算，这条路径上超声波的传播时间为71.488μs。然而，当超声波经过无损伤路径时，传播时间为71.428μs。得到了这种差异之后，可以采用l1最小化方法计算得此处损伤的慢度值为1.0242，识别偏差率为0.8％。最终的识别结果如图18所示，可以保证最大偏差不超过1.7％。需要指出的是，此方法不能完全定位出损伤的精确位置，因为网格的大小大于缺陷的直径，但是可以明确地得出哪个网格内含有缺陷。

情况3

两个裂缝和四个孔洞被设置在真实混凝土结构中，其中孔洞的位置分别为(200，150)，(750，200)，(350，100)和(600，100)，其中两个孔洞的直径为2mm，其另外两个孔洞直径为4mm，稀疏的分布在试件中。如图19所示，两个裂缝大小为100mm×0.5mm，且一个是水平裂缝一个是竖直裂缝。基于惠更斯-菲涅尔原理，超声波绕过孔洞时的路径长度和传播时间计算与上文相同，对于超声波绕过裂缝的情况我们在此详细讨论。如图19所示的竖向裂缝，取不同角度绕行的传播路径，对于1路径，传播距离为300.10mm，2路径的传播距离则为430.52mm。基于传播时间，通过最小化算法可以完成损伤识别，最大误差率为2.3％。需要指出的是，实际路径的传播时间必须和参考的传播时间相比较才能识别出缺陷的位置和形状，这也是本文提出的方法的核心。根据上文提到的裂缝缺陷的各向异性，我们甚至可以不测量最大传播时间和最小传播时间而得到裂缝的方向。

Claims (5)

1.一种基于压缩采样理论的超声CT检测方法，其特征在于，步骤如下：

第一步：确定测量精度，将要监测的试件在监测平面划分成数个单元；

第二步：用计算机生成一个随机的测量矩阵Φ，用于确定测量路线，并用Matlab计算出由各测量路线穿过试件所经过的各单元长度构成字典矩阵A；

第三步：根据上步得到的测量矩阵，从传统的超声CT方法的全测量路径中随机地挑选出本次测量路径；

第四步：按照第三步得出的测量路径在试件两侧放置发射器和接收器，得到测量路径的传播时间

第五步：重复第四步，得到所有选择路径的传播时间，进而构成此次检测的测量向量T′_m；

第六步：建立试件重构模型如下式：

ΔT′_m＝T′_m-T₀＝Φ(A·S-A·S₀)＝Φ·A·ΔS    (1)

其中，ΔT′_m为实际检测的传播时间向量T′_m与相应路径构成的无损伤状态下的标准传播时间向量T的差向量，被检测试件的慢度向量s，其与无损伤状态下的标准慢度向量的差为慢度差向量ΔS，由于损伤在试件中是稀疏分布的，因此该慢度向量ΔS也是稀疏的向量；

第七步：采用基于压缩采样的凸优化方法进行求解，表达式如下

$\mathrm{\&Delta;}\hat{S}=\arg \min ({\left|\right|\mathrm{\&Theta;}\&CenterDot;\mathrm{\&Delta;S}-\mathrm{\&Delta;T}\left|\right|}_2+\mathrm{\&lambda;}{\left|\right|\mathrm{\&Delta;S}\left|\right|}_1)---\left(2\right)$

进而有

$\hat{S}=S_0+\mathrm{\&Delta;}\hat{S}---\left(3\right)$

根据重构出的慢度向量还原试件内部损伤情况。

2.根据权利要求1所述的一种基于压缩采样理论的超声CT检测方法，其特征在于，测量的总次数为m次，而传统方法的测量次数为N次，m＜＜N；随机的测量路径的选择是由测量矩阵Φ来决定的，随机测量矩阵Φ中的所有元素都为0或1，并且每一行只有一个元素为1其余元素均为0，如果第r行第j列的元素为1，意味着选取第j条路径作为第r次随机观测的观测路径，根据测量矩阵Φ从传统的N条全测量路线中随机选取取m条测量路线，其中m满足压缩采样的关系m＞μ·K·log(N/K)，μ因问题不同取值不同，μ≈4，这样试件的内部情况就能够被精确重构。

3.根据权利要求1所述的一种基于压缩采样理论的超声CT检测方法，其特征在于，在测量中，每个路径读三个传播时间，取平均值作为该路径信号的传播时间。

4.根据权利要求1所述的一种基于压缩采样理论的超声CT检测方法，其特征在于，当超声波传播的路径中有孔洞或者缺陷存在时，首波的到达时间会比正常的时间要长，而时间改变量的大小与缺陷的大小成正比，根据不同测量路径的传播时间的改变量，能够确定缺陷的位置和大小。

5.根据权利要求1所述的一种基于压缩采样理论的超声CT检测方法，其特征在于，将监测区划分为p×q＝w个规则网格，w是网格的总数，其中p、q为划分的矩形单元网格数，监测区共有N条射线通过，L_i为第i条射线，T_i是第i条射线的旅行时间，定义速度的倒数为慢度s＝1/V，由Radon公式：

$T_i={\&Integral;}_{L_i}\frac{1}{V_j(x,y)}\mathrm{dL}={\&Integral;}_{L_i}{s}_j(x,y)\mathrm{dL}---\left(4\right)$

其中：

t_ij——是第i条射线在第j个单元中传播的时间

V_j(x，y)——第j个单元中波的传播速度

s_j(x，y)——第j个单元中波的传播的慢度

L_i——第i条射线的长度

设网格足够细，成像单元足够小，每一个单元中的速度V_j(x，y)和慢度s_j(x，y)均视为为一个常数，将上式的积分写成如下累计求和的形式：

$T_i=\underset{j=1}{\overset{w}{\mathrm{\&Sigma;}}}{t}_{\mathrm{ij}}=\underset{j=1}{\overset{w}{\mathrm{\&Sigma;}}}\frac{a_{\mathrm{ij}}}{V_j}=\underset{j=1}{\overset{w}{\mathrm{\&Sigma;}}}{a}_{\mathrm{ij}}{s}_j---\left(5\right)$

探测中，总共发出的N条超声波射线，每一条射线到达对面的行进时间构成了一个时间向量T＝(T₁ t₂ … t_N)^T，形成成如下方程组的形式：

$\left\{\begin{array}{c}{T}_1=\underset{j=1}{\overset{w}{\mathrm{\&Sigma;}}}\frac{a_{1j}}{V_j}=\frac{a_{11}}{V_1}+\frac{a_{12}}{V_2}+...+\frac{a_{1w}}{V_w}\\ T_2=\underset{j=1}{\overset{w}{\mathrm{\&Sigma;}}}\frac{a_{2j}}{V_{j\&CenterDot;}}=\frac{a_{21}}{V_1}+\frac{a_{22}}{V_2^{\&CenterDot;}}+...+\frac{a_{2w}}{V_w}\\ \begin{array}{c}.\\ .\\ .\end{array}\\ T_n=\underset{j=1}{\overset{w}{\mathrm{\&Sigma;}}}\frac{a_{\mathrm{nj}}}{V_j}=\frac{a_{n1}}{V_1}+\frac{a_{n2}}{V_2}+...+\frac{a_{\mathrm{nw}}}{V_w}\end{array}\right.---\left(6\right)$

改写成慢度的形式：

$\left\{\begin{array}{c}{T}_1=\underset{j=1}{\overset{w}{\mathrm{\&Sigma;}}}{a}_{1j}{s}_j=a_{11}{s}_1+a_{12}{s}_1+...+a_{1w}{s}_w\\ T_2=\underset{j=1}{\overset{w}{\mathrm{\&Sigma;}}}{a}_{2j}{s}_j=a_{21}{s}_1+a_{22}{s}_1+...+a_{2w}{s}_w\\ .\\ .\\ .\\ T_n=\underset{j=1}{\overset{w}{\mathrm{\&Sigma;}}}{a}_{\mathrm{nj}}{s}_j=a_{n1}{s}_1+a_{n2}{s}_1+...+a_{\mathrm{nw}}{s}_w\end{array}\right.---\left(7\right)$

改写成矩阵的形式：

T＝AS    (8)

其中：

T——是所有N次检测超声传播的时间向量；

A——是所有N次检测在w个单元中的长度矩阵；

a_ij——超声波在第i个路径上经过第j个单元的长度；

S——所有w个单元的慢度向量。