基于L曲线法的超声CT正则化参数优化

张 培，吕晶晶，李 媛

（中北大学仪器科学与动态测试教育部重点实验室，山西太原 030051）

通过矩量法将上述两个卷积方程转化为代数方程组的形式:

方程（3）、（4）仅描述了物体在某一方向入射波作用下接受到的散射波，为获得足够多的信息，需要从K个方向发射超声波，由此得到K×M个方程，为使上述方程组可解，通

0 引言

80年代初，Johnson将微波技术中矩量法引入超声成像领域［1］，矩量法不再局限于Born近似的限制，可以方便地引入被测物体的某些先验知识（这些先验知识是克服任何逆散射问题的不适定性和对噪声敏感性的关键），可以对大物体、高对比度、强散射的物体进行成像，因此适应范围更广［2］。

不适定性是图像重建逆问题求解中非常普遍的情况，本文采用Tikhonov和TSVD正则化方法解决其不适定性，两种正则化算法的成功应用依赖于正则化参数λ的合适选取。其中容差原理法确定正则化参数需要知道数据项中噪声大小，而广义交叉验证法操作较为繁琐且结果不稳定。本文基于L-曲线法，提出以解的范数和残差变化量的加权形式作为确定正则化参数的依据，在迭代过程根据问题不适定性程度，自适应地调整搜索范围。通过两种正则化技术和优化策略来实现在空间域内的图像重构。

1 UCT图像重建原理

当用超声波从某一方向作用于物体的被成像截面时，超声波与物体内部介质相互作用产生的全场满足非齐次亥姆霍兹方程［3］，其解可用Lippmann-Schwinger积分方程表示:

式中:p（→）和pin分别为全场和入射场;G（→-） 为自由空间的格林函数;O（=k20［n2（→）-1］为物体内部介质声学特性的未知函数。

在实际应用中，我们只能测到物体的散射场和入射场，将式（1）转化为散射场形式:常KM＞N，即要求方程组为超定的。

通过矩量法将上述两个卷积方程转化为代数方程组的形式:

方程（3）、（4）仅描述了物体在某一方向入射波作用下接受到的散射波，为获得足够多的信息，需要从K个方向发射超声波，由此得到K×M个方程，为使上述方程组可解，通

这个超定方程组的求解具有较强的非线性。这里采用Born迭代方法［4］。具体描述如下:

①先假设全场等于入射场P（t）=P（in），由方程（4）求得未知函数的初始解O0，常称其为Born逆解;

②然后将Ok代入方程（3）求更接近实际的物体内部的全场P（kt），k表示迭代次数;

③由第②步求得的全场P（kt），代入方程（4），计算散射场P（ks）=DOP（kt），并求计算散射场与测量散射场的差ΔP（ks）=P（ks）-P（s），若‖ΔP（ks）‖小于事先给定的δ，则停止迭代，否则，转④;

④根据散射场的改变量 ΔP（ks），由方程 ΔP（ks）=DΔOkP（kt求未知函数的改变量ΔOk，然后赋Ok+1=Ok+ΔOk作为未知函数的新值;

⑤由新的未知函数Ok+1（rn），返回步骤②，进一步求更近似的全场。

在迭代过程中，将遇到不适定性逆散射方程的求解问题，而对不适定性问题的正则化技术的选择是否得当，将直接影响迭代算法的稳定性。所以，正则化技术是图像重建问题的关键。

2 正则化方法及参数的选取

算法1:用Tikhonov方法计算Ax=b的算法如下［6］。

①矩阵A的奇异值分解:［U，S，V］=svd（A），其中S=diag（σ1，σ2，…σn）

整数λ称为正则化参数。

③作图（lg‖xλ‖2），lg‖Axλ-b‖2）），确定曲线拐角点，拐点对应的xλ是正则解xtik。

算法2:用TSVD方法计算Ax=b的算法如下［7］。

①矩阵 A 的奇异值分解:［U，S，V］=svd（A）。

②对1≤k≤T，计算解的范数:

③作图（lg（‖xλ‖2），lg（‖Axλ-b‖2）），确定曲线拐角点，拐点对应的xk是正则解xtsvd。

研究表明，对离散不适定问题，lg（‖xk‖2），lg（‖Axk-b‖2）的关系曲线总是呈现L-曲线形状，该关系曲线有一个明显的拐角，最优正则化参数位于L-曲线拐点附近［8］。因此，可通过寻找L-曲线拐点作为最优正则化参数。

3 重建结果比较

实验装置我们采用圆环形结构，发射器（同时也是接收器）等间隔地放在围绕物体的圆环上，半径取20×λ，背景介质为水。采用F=200 kHz的超声波作为入射波，超声波在水中的波速约为C0=1 500 m/s，波长约为7.5 mm，对应的波数为k0=0.837 758 mm-1。按照矩量法，我们对包含物体的正方形区域进行均匀采样，采样间隔均为d=λ/10=0.75 mm，采样点的个数为35×35，划分后单元的总数为N=1 225，且在每个小单元上做内接圆，其半径为a=d/2=0.375 mm。换能器共有50个，可工作于发射、接收两种模式，每次仅有一个换能器发射，其他换能器接收来自物体的散射波。

如图1所示，实验选用的样品如下:对比度为20%的人工绘制的图形，轮廓信息较为丰富。对于已知像函数O（→r），先由全场方程（3）求出ROI中的全场分布，然后通过散射场方程（4）求出物体外的散射场P（s）分布，再使用BI进行图像重建。两种方法的重建结果如图2和图3所示:

图1 目标原始函数

图2 采用Tikhonov迭代4次结果

图3 TSVD正则化方法迭代18次的结果

图4 Tikhonov方法相对误差曲线

由图4可以看出采用Tikhonov方法时，相对误差在第4次相对误差最小，当超过第4次后，相对误差变大。图5可以看出采用TSVD方法迭代到第18次时最小误差恒定不变，因此在第18次中止迭代。对比两种方法，TSVD方法明显好于的Tikhonov方法重建的结果。

图5 TSVD方法相对误差曲线

4 结论

本文运用Born迭代算法解决了方程的非线性问题，对于离散不适定问题，采用Tikhonov和TSVD两种正则化方法分别对方程进行了求解。经实验验证:TSVD正则化方法明显优于Tikhonov正则化方法。超声逆散射成像技术是一个比较复杂的系统工程，算法的稳定性及速度只是问题的一个方面。目前尚没有成熟的理论和技术可用于产品的自动检测，对该方面的一些关键技术进行研究具有一定的学术价值和广阔的应用前景。

［1］Tracy M L，Johnson S A.Inverse Scattering Solutions by a Sinc Basis，Multiple Source，Moment Method-Part II:Numerical Evaluations［J］.Ultrasonic Imaging，1983，5（4）:376-392.

［2］陶进绪，张东文，叶寒生.频域法超声逆散射成像［J］.信号处理，2009，25（2）:169-173.

［3］刘超，汪元美.超声层析成像的理论与实现［D］.浙江大学，生物医学工程，博士论文，2003.

［4］Chew W C，Wang Y M.Reconstruction of Two-Dimensional Permittivity Distribution Using the Distorted Born Iterative Method［J］.IEEE Transactions on Medical Imaging，1990，9（2）:218-225.

［5］韩旭，刘杰，李伟杰.时域内多源动态载荷的一种计算反求技术［J］.力学学报，2009，41（4）:595-601.

［6］Hansen P C，O＇Leary D P.The Use of the L-Curve in the Regularization of Discrete Ill-Posed Problems［J］.SIAM Journal on Scientific Computing，1993，14（6）:1487-1503.

［7］Hansen P C，Jensen T K.An Adaptive Pruning Algorithm for the Discrete L-curve Criterion［J］.Journal of Computational and Applied Mathematics，2007，198（2）:483-492.

［8］王化祥，何永勃，朱学明.基于L曲线法的电容层析成像正则化参数优化［J］.天津大学学报，2006，39（3）:306-309.