摘要_基于射线理论的超声波层析成像技术涉及两个问题_...

黄政宇
摘要：基于射线理论的超声波层析成像技术涉及两个问题：射线追踪技术和方程组的求解技术。

文中对此分别进行了探索：（1）在LTI 射线追踪算法的基础上,采用扩张收缩扫描算法计算最小走时,记录次级源追踪射线路径。

（2）加权阻尼最小二乘反演算法引入充分体现混凝土超声成像特点的权矩阵，以提高反演算法的抗噪能力，使得反演图像唯一、稳定。

对6个模型试件进行速度层析成像的结果表明，本文采用的反演方法可以较好地改善图像质量，有效地反映出混凝土模型试件中缺陷的位置、大小和性质。

最后探讨了影响缺陷成像质量的一些因素。

近年来，在结构混凝土质量检测中引入的超声波层析成像技术，采用图像的方式反映混凝土的内部质量，结果直观明晰，较传统方法有明显的改进。

层析成像（Computerized Tomography ），是在不损伤研究“对象”内部结构的条件下，利用某种场源，根据从“对象”外部用检测设备所获得的投影数据，依照一定的物理和数学关系，利用计算机反演“对象”内部未知的某种物理量的分布，重现“对象”内部特征。

依据不同的原则，层析成像有不同的分类方式：按照射线源类
型，层析成像可分为声波CT 、
射线CT 和电磁波CT 等；根据反演参数的物理性质，可分为波速CT 、衰减CT 和电阻抗CT 等；根据反演的理论基础，可分为以射线理论为基础的图像重建技术和以波动理论为基础的层析方法。

而混凝土无损检测中应用较多的是基于射线理论的超声波速度反演。

本文的结构组织如下：第一节简要阐述基于射线理论的速度反演理论；第二节改进了线性插值射线追踪技术；第三节讨论了求解方程组的加权阻尼最小二乘算法；第四节进行了6个模型试件的速度层析成像，第五节分析了影响缺陷成像质量的若干因素；最后一节给出全文的总结。

1 1 速度反演的原理速度反演的原理速度反演的原理
假定()y x v ,是二维速度模型，如图1所示，v s 1=为慢度，则第i 条射线的初至走时，即投影值，可由式（1）计算
()dl y x s t i
R i ∫=,()N i ,,2,1Λ= (1)
式中i R 是第i 条射线的轨迹，N 为射线条数。

这是一个非线性问题，一般将反演区域离散成若干规则的网格单元（像元），在离散的每个像元内把射线路径近似地看为直线，设ij r 为第i 条射线穿过第
j 个单元的距离，j s 为第j 个单元
的慢度，M 为网格单元总数，则有 ∑==M i j
ij i s r t 1 ()N i ,,2,1Λ= (2)
线性方程组式（2）写成矩阵形式为
T RS = (3)
其中R =(r ij )N ×M 为射线路径矩阵；S =(s j )M ×1为慢度向量；T =(t i )N ×1为走时向量。

求解式（1）为正演过程，即对射线进行追踪，计算走时向量；求解式（3）为反演过程，即获得各像元的慢度s（或速度s v
1=）。

反演主要涉及两个问题，其一为射线路径矩阵的确定，
其二为方程组的求解。

2线性插值射线追踪的改进算法线性插值射线追踪的改进算法
最简单的射线追踪是直射线追踪，即把介质视为均匀完好的，用直线把发射点和接收点连起来，作为射线传播的路径。

当混凝土中存在声速异常区域时，超声波的传播路径会发生明显的弯曲，故宜采用弯曲射线追踪技术。

目前应用效果好的弯曲射线追踪方法有两类：一类是基于程函方程的有限差分法，另一类是基于惠更斯原理的最短路径法。

线性插值射线追踪法(LTI)[5]是程函方程有限差分射线追踪的高级形式。

本文在LTI 的基础上,采用扩张收缩扫描算法,增加次级源的记录，使之适应复杂的速度模型，同时保留了LTI 所具有快速、全局性、编程简单的优点。

算法分两步进行，第一步向前计算最小走时，第二步向后追踪射线路径，其步骤如下： 第一步(向前处理)：计算离散网格模型上所有边界点的最小走时并记录次级源.
a.找到发射点所在单元边界上各计算点的最小走时和次级源；
b.找到发射点单元所在列各单元边界上所有计算点的最小走时和次级源；
c.扩张扫描:扫描由发射点单元所在列开始,向右(左)逐列进行,直到模型的右(左)边界;然后从发射点单元所在行开始向上(下)逐行扫描,直到模型的上(下)边界.扩张扫描过程中,更新最小走时及次级源；
d.收缩扫描:从模型上下!左右边界分别扫描至发射点所在单元的行与列.收缩扫描过程中,
修正最小走时与次级源；
e.重复步骤(c)，再进行一次扩张扫描,则可得每一网格边界点的精确初至走时和相应的次级源。

第二步(向后处理)：追踪所有发射-接收点对的射线路径。

a.首先找到接收点所在单元边界上走时最小的离散点；
b.根据第一步得到的离散点，依据走时的最大梯度,找出其相邻线段上的射线交点；
c.如果所求交点为边界离散点,则找出该离散点的次级源作为新的接收点,否则将射线交点视为新的接收点,重复步骤(a)(b)，直到发射点单元为止；
d.将发射点与最后的射线交点相连,即完成射线路径追踪。

本算法中局部走时的线性插值算法同LTI[6],在此不再赘述.
3 3 二次优化问题及其加权矩阵二次优化问题及其加权矩阵二次优化问题及其加权矩阵
可见，波速反演最终归结为求解大型线性方程组式（3），其基本特点是：射线路径矩阵R 大型、稀疏；方程不相容；条件数高，即走时向量T 和射线路径矩阵R 的误差均会引起解的严重变异。

故方程组式（3）无精确解。

现在为了求得满意解通常要采用某种最优化准则，按l 2范数处理为二次优化问题，并以如下两项准则导出目标函数：
1）最小距离准则：寻求解慢度向量S ，使模2RS T −极小。

2）先验约束准则：对慢度向量附加某种先验约束，以消除满足准则（1）的解的不唯一性。

按上述两准则形成如下二次优化增广目标函数:
()()()()()0201S S W S S RS T W RS T S F T
T −−+−−=µ (4) 式中，1W 为数据加权矩阵，2W 为模型加权矩阵，0S 为先验模型慢度向量，µ为阻尼因子，确定右端两项的相对重要性。

min ()S F
的一阶必要条件为
()()0121RS T W R S W R W R T T −=∆+µ (5) 此即目标函数()S F 的正态方程，其中0S S S −=∆。

方程组式（3）的解为
()()011210RS T W R W R W R S S T T −++=−µ (6)
这种反演方法即是加权阻尼最小二乘法[8]，引入权矩阵的目的是为了利用先验信息，提高反
演算法的抗噪能力，并使得反演图像唯一、稳定。

显然权矩阵的选取是关键，但因为选取的自由度很大，且不存在普遍适用的权矩阵，所以应该具体情况具体分析。

当I W =1（I 为单位矩阵），0=µ，则求解过程为经典的最小二乘法，当I W =1，I W =2，0≠µ则求解过程为阻尼
最小二乘法。

近年来发展了自然权矩阵[7]
，按照对特征值的分析，有 11−=T W D W =2 (7)
矩阵T 的元素为当前模型的各射线的走时；D 的元素为射线穿越任一单元的长度和与相应单元速度乘积，Berryman [9]
称其为射线分布矩阵。

W 1、W 2均为对角矩阵，对称正定。

这种权矩阵的最大特点是表征了基于以下明确物理原理的先验判断：
1） 超声波作为信息载体，随着传播距离的增加衰减增加，波形信号减弱，噪音增加，故射线路径较短者，信噪比较高，也较接近真实路径。

因此对走时较小射线应加重权，故各条射线的权系数为各自走时的倒数。

2） 一条射线对一单元的穿透可视为一次信息“采访”，穿过该单元的射线越多，则“采访”越多，信息也越多，成像更准确，故对射线较密的成像单元应加重权；另外，成像单元的声速越高，则混凝土质量越好，成像更可靠。

所以成像单元的权系数为所有射线穿越该单元的总长度与该单元声速之乘积。

4模型试验模型试验
4.1.1 试验试件试验试件
本试验的对象是6个内部模拟有不同缺陷的素混凝土试件，按试件的平面尺寸分成100cm×100cm 和40cm×40cm 两组，所有试件的厚度均为20cm，其中的缺陷均从上至下贯通。

,a b 中的缺陷均为孔洞，f 中的缺陷为麻石块，g 中为砂浆块，h 中为蜂窝块，i 中为孔洞，如图2所示。

a、b 二试件共给出了5个不同直径的孔洞缺陷，用来研究超声波CT 对缺陷大小和位置的反映情况。

i h g f ,,,四个试件中均只有一个模拟缺陷，用来考察超声波CT 对缺陷性质的反映情况。

初至走时的采集
4.2初至走时的采集
试件平面（即反演区域）进行网格划分，网格（像元）大小为5cm×5cm，发射换能器和接收换能器位于每个网格的中心，通过凡士林作为耦合剂与打磨光滑试件表面紧密接触，移动接收换能器和发射换能器的位置进行超声测试，完成初至走时的采集工作。

测试方式采用双向检测方式，如图3所示。

图中S为发射点，R为接收点。

反演图像
4.3反演图像
利用第二节的射线追踪技术和第三节的加权阻尼最小二乘反演算法进行初至走时反演，然后对反演的原始图像采用方差截断、中值滤波和聚类分析[10][10]
的方法进行图像后处理，得到最终的反演结果如图4所示。

由此可见，本文采用的反演方法可以较好地改善图像质量，有效地反映出混凝土模型试件中缺陷的位置、大小和性质。

5影响缺陷成像的因素分析影响缺陷成像的因素分析
5.1 5.1 射线追踪方式射线追踪方式射线追踪方式
图5 (a)、(b)分别为试件i 直射线追踪和弯曲射线追踪反演的原始图像。

由图可知：直射线追踪和弯曲射线追踪的反演图像类似，两者都能正确反映孔洞缺陷的位置和大小。

直射线追踪反演的混凝土背景声速为4 483m/s，与真实背景声速4 530m/s 的相对误差为-1.04%，孔洞缺陷的声速反演值为3 505 m/s，与背景混凝土声速的相对差异为-22.63%；而弯曲射线追踪反演的混凝土背景声速为4 537m/s，与真实声速的相对误差仅为0.15%，缺陷的声速反演值为3 026m/s，与背景混凝土声速的相对差异为-33.20%。

弯曲的射线路径符合超声波在此试件中的传播特点，因此基于弯曲射线追踪的反演结果有利于反映缺陷的严重程度。

5. 2 5. 2 反演算反演算反演算法法
对于试件b ，当采用最小二乘算法（LS ）和阻尼最小二乘算法（LM ）时，均反演出某些像元的声速为负值，这与实际情况不符，且迭代均发散，无法得到正确的结果。

由于反演均告失败，故未给出最终的层析图像；但当采用基于自然权的加权阻尼最小二乘算法（WLM ），迭代过程没有出现声速参量反演值为负的异常现象，并且计算收敛，反演结果见图4(b)。

5. 3 5. 3 检测方式检测方式检测方式
图6(a)、(b)是试件b 分别在上下单向和左右单向检测方式下经图像后处理的成像结果。

对比图4(b)，可以看出，双向检测方式下反演结果比上下、左右单向检测方式下的结果能更准确地反映缺陷的位置和大小。

这是因为单向检测只能反映缺陷在某一个方向的情况，且射线路径矩阵的条件数高（此例中为1050
），射线覆盖面积小；双向检测能反映缺陷在两个方向上的情况，对缺陷的反映也更全面，而且射线路径矩阵的条件数大大降低（此例中为102），射线覆盖面积大。

5.4 5.4 缺陷最大尺度缺陷最大尺度
缺陷最大尺度 试件a 内部大、中、小三个圆形孔洞缺陷的面积与试件总面积的比例分别为3.14%、2.01%、0.95%。

由图4(a)可知，反演图像能基本准确地反映试件a 内大、中两个缺陷的位置和大小，但
没有反演出试件a 中的小缺陷。

速度反演难以反映试件a 中的小缺陷，一方面是因为基于射线理论的超声波层析成像存在分辨率的极限，其分辨率的最小尺度为波长尺度。

频率50kHz 的超声波在混凝土中的传播速度为4 500m/s 时，其波长为9.0cm，接近试件a 中小缺陷的最大尺度11cm。

另一方面是因为正演测试数据含有噪音，增加了反演的难度。

5.5 5.5 混凝土龄期混凝土龄期混凝土龄期
分别在混凝土龄期一月和二年时对试件f 进行反演。

结果表明，由于混凝土的声速随着龄期的增加而提高，但麻石的声速不随时间变化，故麻石块在混凝土龄期一月时表现为高速异常区，麻石声速比背景混凝土声速高10.21%，而在混凝土龄期二年时表现为低速异常区，麻石声速比背景混凝土声速低0.88%。

6 6 结结 论
结构混凝土超声检测中引入的层析成像技术，能以图像的方式定量完整地反映层析面上混凝土内部质量，弥补“点”上检测的局限，较传统方法有明显的优势，是一种有独特效果的无损检测手段。

由于超声波在混凝土中的传播特性以及混凝土本身的不均匀性，特别是当混凝土中存在蜂窝、孔洞等低速异常体，以及岩石等高速异常体时，超声波会产生折射、绕射等使其传播路径发生明显的弯曲。

因此适合混凝土超声CT 的射线追踪应当是弯曲射线追踪。

加权阻尼最小二乘反演算法引入权矩阵的目的是为了利用先验信息，使反演图像唯一、稳定。

显然权矩阵的选取是关键，但由于不存在普遍适用的权矩阵，所以应该采用体现混凝土超声成像特点的权矩阵，如自然权矩阵。

混凝土中缺陷成像质量还与测试方式、缺陷最大尺度和混凝土龄期相关。