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

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