反射型超声衍射CT迭代重建算法

的正则化目标函数：
ｍｉｎ｛｜｜!"ｆ－ "Ｆ｜｜＋$２｜｜"ｆ｜｜｝
（ ７）
!
其中 $ 为正则化参数， ! 为系数矩阵， 在超声衍射成
像中相应于非均匀傅立叶变换算子， "ｆ＝ ｆｋ＋１′－ ｆｋ′， ｆ０′ 为最初假设原始图像； "Ｆ＝Ｆ′－ Ｆ， Ｆ′是对 ｆ′进行非均
匀傅立叶变换得到的矩阵； 根据傅立叶衍射定理， 由
ｒｉｅｒ ｔｒａｎｓｆｏｒｍ， ｒｅｇｕｌａｒｉｚａｔｉｏｎ
１引言
超声衍射成像是声学成像中的一个重要研究领 域， 除了在医学诊断中的潜在应用外， 还广泛应用于 海 洋 勘 探 、地 震 预 报 、国 防 军 事 、工 业 无 损 检 测 等 领 域。傅立叶衍射定理（ Ｆｏｕｒｉｅｒ ｄｉｆｆｒａｃｔｉｏｎ ｔｈｅｏｒｅｍ） 是 超声衍射成像的理论基础［１］， 该 方 法 是 在 一 阶 Ｂｏｒｎ 近似或者 Ｒｙｔｏｖ 近似前提下， 利用快速傅立叶变换 （ ＦＦＴ） 和反变换（ ＩＦＦＴ） ， 完成层析成像。傅立叶衍射
等。从实际应用角度来看， Ｔｉｋｈｏｎｏｖ 正则化法和截断 奇异值法因为易于实现而应用广泛。
如前所述， Ｔｉｋｈｏｎｏｖ 正则化是通过在正 则 化 目
标函数中增加一个限制条件， 参见式（ ６） ， 适当地选
择限制条件及其在目标函数中的权重， 得到适定解。
本文研究的是反射型超声衍射成像， 由傅立叶衍射
ｎ１＝－ Ｎ１／２ ｎ２＝－ Ｎ２／２
式中（ !１， !２） 代表非均匀频率点的位置。上式用矩阵
表示为：
Ｆ＝"ｆ
（ ２）
Ｙ 是非均匀傅立叶算子。通过 Ｙ 的设定， 可以
得到任意频率点上的傅立叶变换值， 即得到目标物
体二维傅立叶变换域中对应的圆弧上的值。
本文的图像重建， 是利用声场测量数据做一维
傅立叶变换， 得到二维傅立叶变换域上对应的圆弧
图 １ 单频投影频率域内数据点分布图 Ｆｉｇ．１ Ｔｈｅ ｓｉｎｇｌｅ ｆｒｅｑｕｅｎｃｙ ｄｉｓｔｒｉｂｕｔｉｏｎ ｉｎ ｆｒｅｑｕｅｎｃｙ ｄｏｍａｉｎ
２ 算法原理
二维图像的非均匀傅立叶变换定义为［２］
Ｎ１／２－ １ Ｎ２／２－ １
! ! Ｆ（ !１， !２） ＝
ｆ（ ｎ１， ｎ２） ｅｘｐ（ － ｉｎ１!１） ｅｘｐ（ － ｉｎ２!２） （ １）
定理可知， 在频域中得到的数据为高频频率值， 即中
心区域的低频信息几乎没有。Ｔｉｋｈｏｎｏｖ 正则化过程
中保留了所有奇异值， 合理地选择正则化因子， 调节
小的奇异值对重建过程的影响， 可以较好地重建图
像。这也我们是选择 Ｔｉｋｈｏｎｏｖ 正则化方法的原因。
利用 Ｔｉｋｈｏｎｏｖ 正则化进行处理， 可以采用下面
第 ２５ 卷第 ４ 期 ２００６ 年 ０８ 月
声学技术 Ｔｅｃｈｎｉｃａｌ Ａｃｏｕｓｔｉｃｓ
Ｖｏｌ．２５， Ｎｏ．４ Ａｕｇ．， ２００６
பைடு நூலகம்
反射型超声衍射 ＣＴ 迭代重建算法
刘 凌， 冯玉田， 王朔中
（ 上海大学通信与信息工程学院， 上海 ２０００７２）
摘要： 本文以傅立叶衍射定理为基础， 将非均匀傅立叶变换和迭代法相结合， 用正则化方法处理迭代的收敛问题， 建
增加求解过程的稳定性， 用于平衡 #（ Ｆｍ－ !ｆ） 和 Ｌ（ ｆ） 在目标函数中的权重。
（ ３） 采用适当的数值计算方法有效地求解目标 函数， 保证计算的准确性和稳定性。
正 则 化 方 法 很 多 ， 如 Ｔｉｋｈｏｎｏｖ 正 则 化 法 、截 断 奇异值法（ ＴＳＶＤ） 、全变差最小正则化重建法、最大 后 验 正 则 化 重 建 法 、方 差 一 致 限 定 正 则 化 重 建 法 等
接的关系：
Ｆ＝!ｆ
（ ５）
（ ２） 定义正则化目标函数。例如定义一个最小
化目标函数：
"（ ｆ） ＝#（ Ｆｍ－ !ｆ） ＋$Ｌ（ ｆ）
（ ６）
当 "（ ｆ） 取最小值的时候得到 ｆ 的最佳估计。其中 #
（ Ｆｍ－ !ｆ） 用来衡量测量数据 Ｆｍ 和理论结果 !ｆ 之间的 误差， 即重建结果的准确性。引入 $Ｌ（ ｆ） 的目的在于
立了反射型超声衍射成像算法。数据直接在频域中的非均匀频率点上比较， 避免了频域内插引入的误差。本算法也减
少了采样数据量， 降低了运算的复杂度。实验结果表明， 在迭代次数不多情况下， 重建图像可以达到较好的视觉效果。
关键词： 反射超声层析成像； 傅立叶衍射定理； 非均匀傅立叶变换； 正则化
中图分类号： ＴＮ９１１．７
文献标识码： Ａ
文章编号： １０００－３６３０（ ２００６） －０４－０３２６－０５
Ｉｔｅｒ ａｔｉｖｅ ｒ ｅｃｏｎｓｔｒ ｕｃｔｉｏｎ ａｌｇｏｒ ｉｔｈｍ ｆｏｒ ｕｌｔｒ ａｓｏｎｉｃ ｒ ｅｆｌｅｃｔｉｏｎ ｍｏｄｅ ｄｉｆｆｒ ａｃｔｉｏｎ ｔｏｍｏｇｒ ａｐｈｙ
ＬＩＵ Ｌｉｎｇ， ＦＥＮＧ Ｙｕ－ｔｉａｎ， ＷＡＮＧ Ｓｈｕｏ－ｚｈｏｎｇ
（ Ｓｃｈｏｏｌ ｏｆ Ｃｏｍｍｕｎｉｃａｔｉｏｎ ａｎｄ Ｉｎｆｏｒｍａｔｉｏｎ Ｅｎｇｉｎｅｅｒｉｎｇ， Ｓｈａｎｇｈａｉ Ｕｎｉｖｅｒｓｉｔｙ， Ｓｈａｎｇｈａｉ ２０００７２， Ｃｈｉｎａ）
Ａｂｓｔｒ ａｃｔ： Ｂａｓｅｄ ｏｎ ｔｈｅ Ｆｏｕｒｉｅｒ ｄｉｆｆｒａｃｔｉｏｎ ｔｈｅｏｒｅｍ， ａ ｒｅｆｌｅｃｔｉｏｎ ｍｏｄｅ ｄｉｆｆｒａｃｔｉｏｎ ｔｏｍｏｇｒａｐｈｙ （ ＲＭＤＴ） ｒｅｃｏｎｓｔｒｕｃｔｉｏｎ ａｌｇｏｒｉｔｈｍ ｉｓ ｅｓｔａｂｌｉｓｈｅｄ ｂａｓｅｄ ｏｎ ｔｈｅ ｃｏｍｂｉｎａｔｉｏｎ ｏｆ ｔｈｅ ｎｏｎ－ｕｎｉｆｏｒｍ ｆａｓｔ Ｆｏｕｒｉｅｒ ｔｒａｎｓｆｏｒｍ （ ＮＵＦＦＴ） ａｎｄ ａｎ ｉｔｅｒａｔｉｖｅ ｍｅｔｈｏｄ． Ａ ｒｅｇｕｌａｒｉｚａｔｉｏｎ ｍｅｔｈｏｄ ｉｓ ｕｓｅｄ ｔｏ ｅｎｓｕｒｅ ｃｏｎｖｅｒｇｅｎｃｅ ｏｆ ｔｈｅ ｉｔｅｒａｔｉｏｎ． Ｗｉｔｈ ｄｉｒｅｃｔ ｃｏｍｐａｒｉｓｏｎ ｏｆ ｎｏｎ－ｕｎｉｆｏｒｍｌｙ ｄｉｓｔｒｉｂｕｔｅｄ ｆｒｅｑｕｅｎｃｙ ｓａｍｐｌｅｓ ｉｎ ｔｈｅ ｆｒｅｑｕ－ ｅｎｃｙ ｄｏｍａｉｎ， ｅｒｒｏｒｓ ｄｕｅ ｔｏ ｔｈｅ ｆｒｅｑｕｅｎｃｙ ｄｏｍａｉｎ ｉｎｔｅｒｐｏｌａｔｉｏｎ ａｒｅ ｒｅｄｕｃｅｄ． Ａｌｓｏ， ｔｈｅ ａｍｏｕｎｔ ｏｆ ｓａｍｐｌｅｄ ｄａｔａ ｎｅｅｄｅｄ ｉｎ ｔｈｅ ｒｅｃｏｎｓｔｒｕｃｔｉｏｎ ｉｓ ｄｅｃｒｅａｓｅｄ ａｎｄ ｃｏｍｐｕｔａｔｉｏｎ ｃｏｍｐｌｅｘｉｔｙ ｏｆ ｔｈｅ ａｌｇｏｒ－ ｉｔｈｍ ｒｅｄｕｃｅｄ． Ｅｘｐｅｒｉｍｅｎｔａｌ ｒｅｓｕｌｔｓ ｓｈｏｗ ｔｈａｔ， ｗｉｔｈ ａ ｆｅｗ ｉｔｅｒａｔｉｖｅ ｓｔｅｐｓ， ｔｈｅ ｒｅｃｏｎｓｔｒｕｃｔｅｄ ｉｍａｇｅ ｈａｓ ａｃｃｅｐｔａｂｌｅ ｖｉｓｕａｌ ｑｕａｌｉｔｙ． Ｋｅｙ ｗｏｒ ｄｓ： ｒｅｆｌｅｃｔｉｏｎ ｍｏｄｅ ｄｉｆｆｒａｃｔｉｏｎ ｔｏｍｏｇｒａｐｈｙ； Ｆｏｕｒｉｅｒ ｄｉｆｆｒａｃｔｉｏｎ ｔｈｅｏｒｅｍ； ｎｏｎ－ｕｎｉｆｏｒｍ ｆａｓｔ Ｆｏｕ－
ｆ ｋ＋１＝ｆ ｋ＋"ｆ ｋ
（ ９）
每次迭代后计算式（ ７） 的目标函数， 如果大于设定的门 限， 则利用式（ ８） 更新 ｆ， 这样的过程重复进行， 直到求 解的目标函数小于设定的门限值为止， 此时认为迭 代过程收敛， 这时的图像 ｆ ＊ 即被认为是重建图像。
在迭代过程中， 正则化参数 $ 对结果影响很 大， 对于 $ 的选择是正则化的一个难题， 目前关于 正则化参数的选择方法有基于离差原理（ ｄｉｓｃｒｅｐａｎ－ ｃｙ ｐｒｉｎｃｉｐｌｅ） 的方法、广义交叉验证方法（ ｇｅｎｅｒａｌｉｚｅｄ ｃｒｏｓｓ－ｖａｌｉｄａｔｉｏｎ， ＧＣＶ） 、准 最 优 化 （ ｑｕａｓｉ－ｏｐｔｉｍａｌ） 方 法和 Ｌ 曲线（ Ｌ－ｃｕｒｖｅ） 法等， 但具体选择哪种方法选 择参数， 需要根据具体的模型来尝试。其中 Ｌ 曲线 方法［５， ６］ 是比较接近超声逆散射模型， 也是比较常用 的， 因此采用此方法。
上的值， 所以图像重建的根本问题在于求解式（ ２） 的
逆问题， 即
ｆ＝" ＋Ｆ＝（ " Ｈ"） －１" Ｈ"
（ ３）
其中 " ＋为 " 的伪逆。可以看出， 直接求解式（ ３） 的
计算量是很庞大的。如果 " 是一个 Ｎ 个元素的一维
矩阵， 需要 Ｎ×Ｎ 的逆矩阵， 算法的计算复杂度为 Ｏ
（ Ｎ３） ； 如果是二维 Ｎ×Ｎ 矩阵， 需要 Ｎ２×Ｎ２ 的逆矩阵，
克服上述缺点， 不直接求解非均匀傅立叶变换， 对得
到的测量数据不进行插值处理， 而采用正则化迭代
方法重建图像。
如前所述， 原始图像 ｆ 与 Ｆ 之间有式（ ２） 的映射
关系， 实际应用中， 若只知道 Ｆ 的测量数据 Ｆｍ， 就可 以由 Ｆｍ 出发， 得到原始图像 ｆ 的近似估计 ｆ ′， 即重 建图像。因此， 需要找到一个将重建目标的估计 ｆ′
收稿日期： ２００５－１１－２８； 修回日期： ２００６－０２－２７ 作者简介： 刘凌（ １９８１－） ， 女， 浙江人， 硕士研究生， 研究方向： 图形处理。
定理有两种模式， 即透射型超声衍射投影定理和反 射型超声衍射投影定理。透射型衍射是在前向散射 场中采集数据， 它的一维傅立叶变换给出了物体 Ｏ（ ｘ， ｙ） 的二维傅立叶变换 Ｆ（ ｗ１， ｗ２） 的内半侧圆弧上 的值（ 如图 １（ ａ） 所示） 。反射型衍射是在后向散射场中 采集数据， 它的一维傅立叶变换给出了物体 Ｏ（ｘ， ｙ） 的 二维傅立叶变换 Ｆ（ ｗ１， ｗ２） 的外半侧圆弧上的值（ 如 图 １（ ｂ） 所示） 。可见傅立叶衍射投影定理将散射场 采集数据的一维傅立叶变换和被测物体断面的二维 傅立叶变换建立了关系。由于 ＩＦＦＴ 必须由均匀分布 的频率点得到， 所以重建过程包含从非均匀频率点 插值到均匀网格点的过程， 这样就会引入误差。非均