（信号与信息处理专业论文）nufft及频谱外推在超声衍射层析成像中的应用研究.pdf

摘 要 i 摘摘 要要 超声 ct 成像具有无辐射、成像方式更多样、价格便宜等优点在医学成像、 工业无伤检测等领域具有广泛的应用前景。 傅里叶衍射层析成像是建立在傅里叶 衍射定理基础上的成像方法，这种方法基于波动方程，考虑了物体内部介质对入 射声波的散射作用，对物体周围的散射场进行成像，能更全面的反映物体内部信 息， 获得更加准确直观的重建图像。但是建立在傅里叶衍射定理基础上的重建方 法存在一些问题：其一是得到的频域采样数据分布在频域的一些圆弧上，这就引 出了非均匀数据的离散傅里叶变换（ndft）问题，以及如何获得对非均匀数据的 快速傅里叶变换（nufft） ；其二是这些采样点都局限在频率域的低频部分（透射 型超声 ct） ，高频信息的缺失会影响到重建图像的质量。本文主要对这两个问题 做出了分析，给出了解决方法，主要内容如下。 首先针对得到采样数据分布不均匀的问题，本文采用了一种 nufft 的方法。 利用近似理论用一组序列来近似傅里叶基，导出了一种解决 ndft 的一种近似表 达式。 由此获得了一种 nufft 的算法并给出了实现方法。这种 nufft 算法一般包 含插值、过采样的 fft 以及加权处理三个步骤。其中插值函数在整个过程中涉及 到插值和加权处理两个步骤，它的选择直接影响到算法的性能。本文采用了窗函 数 gaussian 窗和 kaiser-bessel 窗来近似，并分别对两种情况进行了误差分析 和计算复杂度的分析，实验表明 kaiser-bessel 窗能得到较优的结果。另外针对 有噪声污染的图像， 介绍了一种利用总分正则化的迭代方法来改进图像质量的方 法。 通过实验验证，在相同条件下使用本文的 nufft 方法能获得比 fbp 及双线性 插值更高质量的重建图像，并且随着采样数的提高能进一步改善图像质量。 针对频率域的采样数据都局限在低频部分的问题。本文首先利用 g-p 算法， 一种利用物体空间受限先验知识的迭代算法， 从 nufft 重建的初始图像谱外推获 得初始的高频信息，再使用非线性谱外推的方法来增强图像的高频部分，得到改 善的图像。另外，本文还对非线性谱外推的方法做出了改进。针对原方法不能使 图像中不同幅度的分量均匀增强的问题， 本文采用了一种能反映图像幅度局部信 息变化的参量来代替原方法中的全局阈值，以使图像能够得到均匀增强。 最后本文将以上算法结合起来， 给出了频域数据的采样开始一直到最后得到 增强图像的重建过程的处理流程。 关键字关键字：超声衍射层析成像 傅里叶衍射定理 nufft g-p 算法 非线性谱外推 abstract ii abstract ultrasound computerized tomography(ct) have many advantages, such as non-invasion, low-cost and so on. and it has been widely used in many applications like medical imaging and industrial non-invasion detection. in the future it will be applied in more and more areas. the fourier diffraction ct imaging method, which uses the fourier diffraction theorem, bases on wave equation. this method takes the scatter field into consideration. the results, reconstructed from the scattered filed of ultrasound wave, carry more information about the object. so the reconstruction results are higher definition and more perceived. however, there are some issues about this reconstruction method. one is that the samples in the frequency domain distribute on some arcs. so they are not uniform distributed. the other problem is that these samples all limited in low frequency band (the transmission type). in this paper, we analyze these two issues and give solutions. first, we solve the issue about the non-uniform distribution samples. we apply a type of nufft method. we apply the approximation theorem to use a sequence to approximate the fourier basis. and we get an approximate representation to solve the ndft. so we now have a type of nufft algorithm and its implementation. the method includes three steps: interpolation , computing over-sampled data using fft， scaling. among these steps, interpolation function is involved in interpolation process and scaling process, so the choice of interpolation function is the key to the performance of algorithm. we select the gaussian window and the kaiser-bessel window as the approximate solution. and we estimate the errors and analyze computation complexity. the simulation shows the kaiser-bessel window can get a more promising result. and for the reconstruction in present of noise, we introduce an iterative reconstruction, making use of total variation regularization to improve the quality of the results. the simulation verifies that under the same conditions the method applied in this paper can achieve a more acceptable result compared with fbp and bilinear interpolation method. second, we discuss the problem that all samples in frequency domain confined in the low-pass part. at first, we make use of g-p algorithm to get high frequency information from the image which is reconstructed using nufft method early discussed. the g-p method is an iterative procedure, which combine information abstract iii about the fourier transform of a function with independent space domain constraints. and then a modified non-linear extrapolation method is applied to enhance the high frequency part. compared with the original method, the improved method replaces the global threshold with a parameter which reflects the local information about amplitudes changes in the image. so the new method can enhance the image uniformly. at last, combining these methods, we build a procedure which can reconstruct a high quality image from projections. keyword: ultrasound diffraction computerized tomography imaging, the fourier diffraction theorem, nufft, g-p algorithm, nonlinear reconstruction. 中国科学技术大学学位论文原创性声明 本人声明所呈交的学位论文,是本人在导师指导下进行研究工作所取得的 成果。除已特别加以标注和致谢的地方外，论文中不包含任何他人已经发表或 撰写过的研究成果。与我一同工作的同志对本研究所做的贡献均已在论文中作 了明确的说明。 作者签名：_ 签字日期：_ 中国科学技术大学学位论文授权使用声明 作为申请学位的条件之一，学位论文著作权拥有者授权中国科学技术大学 拥有学位论文的部分使用权，即：学校有权按有关规定向国家有关部门或机构 送交论文的复印件和电子版， 允许论文被查阅和借阅， 可以将学位论文编入 中 国学位论文全文数据库等有关数据库进行检索，可以采用影印、缩印或扫描 等复制手段保存、汇编学位论文。本人提交的电子文档的内容和纸质论文的内 容相一致。 保密的学位论文在解密后也遵守此规定。 公开 保密（_年） 作者签名：_ 导师签名：_ 签字日期：_ 签字日期：_ 第 1 章 绪论 1 第第 1 章章 绪论绪论 1.1 引言 计算机断层技术是一种通过投影来计算物体的断层切片图像的一种方法。 其 中包括了大量的无损成像技术，大量应用于地质勘探，工业无伤检测等领域，尤 其在医学中以诊断为目的的应用。最早的断层成像的概念可以追溯到 abel,早在 1826 年他就第一次提出从投影重建物体的方法，这种方法利用了投影的几何对 称性。 这种方法在 1917 年由 radon 进行了扩展， 应用到具有特定形状的物体上。 而这些理论的最早的实际应用则诞生于 1972,hounsfield 和 cormack 发明了第一 台 ct 扫描仪，他们也由此分享了 1979 年的诺贝尔奖。之后，发展出了各种类 型的 ct 设备，ct 技术也得到了极大的进步。 记被测物体对 x-射线的线性衰减系数为),(yx，有： dlyx l eii ),( 0 （1.1） 可以得到 )ln(),( 0 i i l dlyx （1.2） ct 图像重建就是根据以上式子各个方向的线积分值（即射线投影）来确定物体 内部组织结构，即反演分布）（ yx,。 x-ct 的诞生以来，获得了极大的发展，主要的方法可以分为变换法和级数 法。其中，常见的变换法有 radon 反变换方法，卷积反投影（fbp）,双线性内 插法等等。卷积反投影算法先对投影函数进行修正，然后在直接进行反投影重建 图像。这种方法精度比较高，但是计算量交大而速度慢。双线性内插算法则直接 运用傅里叶反变换的 快速算法，因而速度较快，但是由于插值的误差的存在导 致成像效果不是很理想。但由于计算速度上的优势，这种方法得到极大的关注， 衍生出很多新类型，具有很大的潜力。级数展开方法主要有代数重建方法和基于 em 算法的极大似然估计重建算法等。代数重建算法（art）最早由 gordan 等 人于 1970 年提出的。代数重建方法是图像重建的问题重新表述为大型方程组的 求解问题。这种方法相比于之前的变化法，重建的图像精度高，但具有运算速度 比较慢的确定，且不易于硬件实现。不过随着计算机技术的发展，这种方法越来 越受到重视。例外一种级数展开法由 lange 等人提出，他们将 em 算法应用到透 射型 ct 中，基本思想是将探测器接收到的光子数目看成是服从泊松分布的随机 变量， 重建问题就变成了获得这样的图像向量使得概率测量数据的团体可能性函 第 1 章 绪论 2 数去的取得最大值。 这种方法特别适用于获得投影数据不完全时重建图像的情况 1.2 超声 ct 及相关技术的发展 超声 ct 由 x-ct 技术发张而来的，相对于 x-ct 而言，超声 ct 具有一定 的优点： （1）安全性。x-ct 的射线源是 x-ray,，具有一定的辐射性，特别在医疗成像应 用中大剂量的使用对会对人体照成一定的损害， 相对来说以超声为检测源的超声 ct 就很安全。 （2）成像方式的多样性。x-ct 只是利用 x-ray 通过物体之后衰减的幅度，然后 对其求对数来作为重建图像的像素值， 超声 ct 可以通过测量不同参量， 如声速、 声衰减系数、剩散射系数以及非线性物理参量来定量的重建图像，不仅能反映生 物组织的形态变化，还能通过这些不同参量直接反映组织的生化特性。 （3）经济性。相对于 x-ct 昂贵的设备，超声 ct 设备价格低很多，同时体积 小，携带方便具有更广阔的应用前景。 也正是超声 ct 具有的以上优点，使之具有广阔的应用前景 超声 ct 由 x-ct 技术衍射而来，早期的方法借鉴了 x-ct 的方法。但是把 超声在物体中的传播路径看成是直线的直线传播的模型， 这对 x-ct 来说是合适 的，相对于物体组织来说 x-ray 的尺寸是很小的，因此一般情况下可以忽略散射 现象。 然而超声的波长较长，当其在不均匀的物体中传播时会有折射及衍射等现 象（统称为散射现象） ，所示将其传播路径看成是直线会造成伪迹。下面简单介 绍下超声成像的一些方法： （1）代数重建方法（art）和射线跟踪方法(ray-tracing approach) 代数重建方法是把切片的折射系数看成是一组未知项， 然后利用测量得到的 投影数据建立方程组，将未知项和投影数据联系起来。这样重建图像的问题就转 变成对方程组的求解问题。虽然概念上相对简单，但是这种方法在实际应用中， 尤其在医疗成像中缺乏准确性， 方程组的解的迭代求解过程可能收敛很慢而导致 速度很慢。不过在某些情况，当我们只能获取部分角度的投影或者投影不是均匀 分布的时候，结合这种方法和其他方法（fbp）就可以帮助我们获取比较理想的 重建图像。 在这种方法中，对于折射系数和投影数据之间的关系，对 x-ct 来说，这个 可以比较简单的可以视作是直线传播的。对于超声 ct 却不是这样。由于物体组 织内部是不均匀的，声波传播路线会发生弯曲。这样我们需要可以先初始的认为 声波在物体组织直线传播的， 然后运用光学中的射线传播方法来修正这一传播路 第 1 章 绪论 3 径， 得到折射系数和投影值之间更加精确的关系式，迭代这个过程直到得到比较 理想的结果。这样，结合以上两种方法可以得到比较精确的重建图像。但是这类 方法一般适用于非均匀小于 23%的物体组织，如果物非均匀度大于 10%，散射 现象将很严重，此类方法得到的结果将失去意义。 （2）衍射层析成像方法 这类方法就是研究弱散射的条件下，充分考虑了介质的不均匀性，及对声波 传播的影响，建立了介质参量和散射声场边界值（接受数据）之间的关系。通过 求解以下非均匀介质条件下的波动方程： )()()()( 2 0 2 ruroruk （1.3） 来得到超声波在非均匀介质中的传播的平面分布。在 born 或者 rytov 近似 下， 可以得到傅里叶衍射定理，建立了投影的一维傅里叶变换和像函数频域之间 的对应关系，由此来重建图像。衍射层析成像的方法充分考虑了散射的影响，具 有较高的精度，同时在计算过程可以引入 fft 算法，极大地提高运算速度，可 以获得实时的图像。因而，这类方法获得极大的关注，发展很快，本文也将主要 讨论此方法。建立在衍射层析成像的方法有基于精确场描述的，通过迭代切结全 场方程和散射方程，来逼近介质像函数以及感兴趣区域（roi）的全场分布。这 种方法本质上来说是将问题看成方程组的求解问题， 因而具有很强的不适定性和 非线性，对于要求分辨率很高的情况，成像的图像尺寸会很大，方程组的规模也 会很大，这样问题复杂度会很高，运算效率就会很变低。 此外超声 ct 还可以被分为透射型和反射性（urct）两种。反射式超声 ct 发射和接受超声波的传感器位于介质的同一侧， 通过接收超声回波信息来重建图 像。这种方法一方面可以克服 b 型超声成像中声线传播方向上排列在后方的物 体在图像中易形成阴影的不足，且能同时得到高的纵向和横向分辨率；另一方面 可以使成像系统简单而易于实现，因而受到人们的重视。研究人员这种方法运用 于超声无损检测，生物组织成像、行星的观察等等应用中。透射型超声发射和接 受超声波的传感器位于介质的俩侧，这也是本文研究的主要方式。 目前，超声 ct 已经在医学成像领域获得极大的应用，并逐步拓展到工业材 料的无损检测，航空航天，军事领域，同时在矿物勘探、地震构造、地震预测预 报等领域也将有着极大的前景。 超声 ct 总的趋势将是像高速，高分辨率、更加安全可靠的方向发展，这就 是数据采集更加方便、快速，成像过程时间更短，等到的重建图像具有更高的分 辨率，更加可靠。此外在不完全投影以及三维图像的重建也是未来超声 ct 的重 大的研究方向。 第 1 章 绪论 4 1.3 本文的研究内容和结构 1.3.1 本文的研究内容和创新 本文主要针对由傅里叶衍射定理重建图像面临的两个问题展开了讨论， 其一 是傅里叶衍射层析成像方法得到的频域采样数据分布在频域的一些圆弧上， 这就 引出了非均匀数据的离散傅里叶变换（ndft）问题，以及如何对非均匀的数据 进行快速傅里叶变换（nufft）的问题；第二个是这些采样点都会局限在频率 域的低频部分（透射型超声 ct） ，高频信息的缺失会影响到重建图像的质量。本 文对这两个问题做出了分析，给出了解决方法，并形成了从傅里叶衍射定理得到 的采样数据到获得最终高质量重建图像的处理过程。 首先针对得到采样数据分布不均匀的问题，本文先对 ndft 的问题进行了 一般性的阐述，接着引入 nufft 的概念。然后利用近似理论，从 shannon 定理 出发，用一组序列来近似傅里叶基，导出了一种解决 ndft 的一种近似表达式。 由此获得了一种 nufft 的算法并给出了实现方法。这种 nufft 算法一般包含 对插值、过采样的 fft 以及加权处理三个步骤。其中由于插值函数在整个过程 中涉及到插值和加权处理两个步骤，它的选择直接影响到算法的性能。本文采用 了窗函数 gaussian 窗和 kaiser-bessel 窗来近似，并分别对两种情况进行了误差 分析和计算复杂度的分析， 实验表明 kaiser-bessel 窗能得到更优的结果。 另外针 对有噪声污染的图像， 介绍了一种利用总分正则化的迭代方法， 来改进迭代质量。 最后进行了仿真实验，并与第二章中介绍的一些基本方法做出了比较。我们发现 在相同条件下使用本文的nufft方法能获得比fbp及双线性插值更高质量的重 建图像，并且随着采样数的提高能进一步改善图像质量。 接着本文对第二个问题， 即频率域的采样数据都会局限在低频部分的问题展 开了讨论。本文首先利用 g-p 算法，一种利用物体空间受限先验知识的迭代算 法，由 nufft 重建的初始图像谱外推获得初始的高频信息，再使用非线性谱外 推的方法来增强图像的高频部分，得到改善的图像。另外，我们这里还对非线性 谱外推的方法做出了改进。 针对原方法不能使图像中不同幅度的分量均匀增强的 问题， 我们这里采用一种能反映图像幅度局部信息变化的参量来代替原方法中的 全局阈值，以使图像能够得到均匀增强。然后我们将整个重建过程结合起来，给 出了频域数据的采样开始一直到最后得到增强图像的处理流程。 1.3.2 本文的结构 在第二章中，我们从非齐次波动方程出发，利用 born 近似以及 rytov 近似 第 1 章 绪论 5 对方程的求解过程进行了简单的推导，最后给出了方程的解析解的表达式。再根 据解的表达式进一步推出了傅里叶衍射定理， 建立投影数据一维傅里叶变换同物 体像函数之间的关系，这将是我们本文所有算法的理论基础。接着我们还讨论了 利用傅里叶衍射定理重建的一些问题。 简单描述了平行投影和扇形投影两种投影 的扫描方式；分析了保证成像质量时，扫描步长需要满足的条件，并对分辨率和 成像质量的关系进行了探讨。 还对成像过程中从理论到实际应用中可能出现的误 差和噪声问题进行了分析。 在第三章中，我们简单介绍了现有的傅里叶衍射层析成像的基本方法。这些 方法大体可以分为两大类。其中一类是频域插值方法，这种方法先利用过采样技 术再通过插值的方法来重建图像，插值函数可以选用近邻插值和双线性插值。另 外一种均匀频域插值方法则将待重建物体的空域受限的性质考虑到插值函数的 推导过程中的。另外一类是空域插值方法，如滤波反投影算法。这两类算法都有 各自的缺点，对于频域插值的方法来说， ，最大的问题在于各种插值方法都会带 来很大的插值误差，造成重建得到的图像不是很精确。但同时有个很大优势在于 可以在计算中利用到快速傅里叶变换（fft） ，从而大大提高重建的速度。而滤 波反投影的方法在能获得足够投影数据的时候可以得到较好的重建结果， 但是其 计算复杂度比较高，而且在不完全投影的情况下，不能得到高质量的重建图像。 后续章节将会与这里方法的进行对比。 在第四章中， 我们针对傅里叶衍射层析成像方法得到的频域采样值分布的不 均匀性问题展开分析。 首先对 ndft 的问题进行了一般性的阐述， 接着引入 nufft 的概念。 然后利用近似理论， 从 shannon 定理出发， 用一组序列来近似傅里叶基， 导出了一种解决 ndft 的一种近似表达式。由此获得了一种 nufft 的算法并给出 了实现方法。这种 nufft 算法一般包含对插值、过采样的 fft 以及加权处理三个 步骤。 其中由于插值函数在整个过程中涉及到插值和加权处理两个步骤，它的选 择直接影响到算法的性能。采用了窗函数 gaussian 窗和 kaiser-bessel 窗来近 似， 并分别对两种情况进行了误差分析和计算复杂度的分析。针对有噪声污染的 图像， 介绍了利用总分正则化的迭代方法来改进迭代质量的方法。最后进行了仿 真实验，并与第二章中介绍的一些基本方法做出了比较。我们发现在相同条件下 使用本文的 nufft 方法能获得比 fbp 及双线性插值更高质量的重建图像， 并且随 着采样数的提高能进一步改善图像质量。针对分别采用 gaussian 窗和 kaiser-bessel窗的两种情况， 仿真结果表明 kaiser-bessel窗具有稍优的性能。 在第五章中， 主要针对建立在傅里叶衍射定理重建方法得到的频域采样数据 局限在低频域的问题展开了讨论。 首先利用 g-p 算法对 nufft 重建的初始图像谱 外推获得初始的高频信息，再使用非线性谱外推的方法来增强图像的高频部分， 第 1 章 绪论 6 得到改善的图像。另外，我们这里还对非线性谱外推的方法做出了改进。针对原 方法不能使图像中不同幅度的分量均匀增强的问题， 采用一种能反映图像幅度局 部信息变化的参量来代替原方法中的全局阈值，以使图像能够得到均匀增强。并 对以上内容进行了仿真实验，实验结果表明改进的算法获得了较优的结果。 第六章对本文进行了总节和展望。 第 2 章 衍射层析成像的基础理论 7 第第 2 章章 衍射层析成像的基础理论衍射层析成像的基础理论 当使用超声 ct 来重建图像时，因为超声波的波长和物体内部介质大小是可 比拟的，所以超声波在物体内部传播过程中的散射现象就不能被忽略。于是就不 可以再简单的应用直线传播的模型了。这就是说，得到的投影数据并不能简单的 表示为直线上的积分，它会受到物体折射系数的影响。某种意义上来说，用几何 声学的来模拟声波传播路径的分析方法将不再适用。 可以采用一些方法来修正这 一问题。其中一种是代数重建方法，首先对折射系数进行一个初始的赋值从而对 声波传播路径进行初始的估计，然后不断迭代来修正折射系数值和传播路径，最 后会收敛到实际的折射系数。这种方法适用于散射非常弱的情况下，如物体介质 的不均匀度小于 3%的时候。但当介质的不均匀程度比较大的时候，这些使用几 何声学理论的方法就会不再适用。 本文这里就会转而使用波动声学的理论来解决 这个问题。本章接下来就从波动方程出发，得到傅里叶衍射定理，进而以此为基 础得到衍射层析成像的方法。最后简单讨论了超声 ct 重建中的一些具体问题。 2.1 非均匀介质中的波动方程及其解 对于声波在物体介质中的传播过程，由于需要考虑散射作用，声波不再沿着 直线传播并且投影值也不再是直线积分，因而不能采用几何声学的模型。所以我 们这里用波动方程来描述这个能量的流动。本节将从理想流体中的波动方程出 发，导出 born 和 rytov 近似条件下的解，最后推出傅里叶衍射定理。 2.1.1 波动方程 在超声 ct 中我们，我们将更关心非均匀介质条件下的波动方程，我们将通 用形式的波动方程表示如下kak 1999： 0)()( 2 rurk （2.1） 这里我们将)( rk视作表示介质折射系数的标量函数，写作： )(1 )()( 00 rnkrnkrk （2.2） 其中， 0 k表示介质的平均波数，)( rn表示折射系数的变化量。一般来说，我们 认为物体限定在有限的范围之内，因此)( rn在物体之外可视为取 0。重写波动 第 2 章 衍射层析成像的基础理论 8 方程如下： )( 1)()( 22 0 2 0 2 rurnkruk （2.3） )( rn表示介质 r处的复数折射系数，等于： )( )( 0 rc c rn （2.4） 这里 0 c代表声波在介质周围物质中的传播速度，)( rc是物体中 r处的声速。 对声 波来说，如果只涉及在可压缩的粘性流体中的情况，我们这里有： )()( 1 )( rr rc （2.5） 和分别代表 r处的局部密度和压缩比。 值得注意的是， （2.3）式是在我们忽略了一阶及高阶有关介质参量微分项的 情况下成立的。 如果我们考虑这些高阶微分项， 波动方程的精确表达式可表示为： )()( 2 0 2 0 2 uukruk k ）（ （2.6） 这里 0 0 k （2.7） 0 （2.8） 0 和 0 根据成像过程的建模方式的可以选择为介质所在环境的压缩比和密度， 或者介质本身的平均压缩比和平均密度。 通过以上的分析，我们可以得到一个更加通用的波动方程的表示形式： )()()( 2 0 2 ruroruk ）（ （2.9） 其中 1)()( 22 0 rnkro （2.10） 我们一般认为)( ru由两部分组成，)( 0 ru和)( rus。其中，)( 0 ru作为入射场，是 不存在非均匀性的情况下的场强，等价于以下齐次方程的解： 0)( 0 2 0 2 ruk ）（ （2.11） )( rus部分作为散射场，作为总场强的一部分可以认为由介质的非均匀引起的。 如上所述)( 0 ru可以作为齐次方程的解，而)()()( 0 rururu s ，将这两部分带 入公式（2.9） ，我们得到只有散射场)( rus的方程： 第 2 章 衍射层析成像的基础理论 9 )()()( 2 0 2 ruroruk s ）（ （2.12） 这种 helmholtz 等式并不能直接求解得到)( rus，但是可以得到由格林函数组成 的解mor53。格林函数是以下微分方程的解： )()|() 2 0 2 rrrrgk（ （2.13） 在三维空间里可以写作： r e rrg rik 4 )|( 0 （2.14） 其中 rrr （2.15） 二维空间里，公式（2.13）的解可以用一阶第一类 hankel 函数表示为： )( 4 )|( 0 )1( 0 rkh j rrg （2.16） 两种情况里面，格林函数是差分 rr 的函数，所以这里我们简单将其表示 为 )( rrg 。因为公式（2.13）表示一个点的不均匀度，所以这里的格林函数也 就表示了一个点的散射场。这样我们有下面的表达式： )()()()()(rdrrruroruro （2.17） 这里我们将不均匀介质中的波动方程表示成不同相移 r的冲击函数同 )()( ruro乘积的和值。格林函数代表单个 delta 函数条件下波动方程的解，由于 波动方程的左边是线性的， 所以我们可以将每个单独的散射点产生的散射场叠加 起来获得方程的解。 通过这种思想，由脉冲)()()( rrruro引发的总场强可以写作脉冲响应 )( rg的权重和位移的和。这是一个简单的卷积过程并且公式（2.12）右边所代表 的所有源点的散射场可以写作一下的叠加： )()()()(rdrurorrgrus （2.18） 这样，物体内部的全场就可以表示为： 00 )()()()()()()(rdrurorrgrurururu s （2.19） 这个表达式称为全场方程，公式（2.18）称为散射场方程，结合这两个方程 式我们研究逆散射问题的基础。通常情况下，当我们已知物体内部结构及其声学 特性参数时， 研究超声波穿过物体时物体的全场及物体外部的散射场的问题称作 前向散射问题。相反的，逆散射问题指的是在已知超声波穿透物体的散射场的情 第 2 章 衍射层析成像的基础理论 10 况下来反演物体内部构造的问题，这也是本文研究的主要问题。接下俩我们讨论 通过两种近似来求解这个方程。 2.1.2 born 近似 born 近似是两种方法中相对简单的一种。 将总场)( ru看成是入射场)( 0 ru和 散射场)( rus的两部分的和，积分(2.18)可以重新表示如下： 0 )()()()()()()(rdrurorrgrdrurorrgru ss （2.20） 这里假如散射场)( rus远小于)( 0 ru， 上式积分中的第二项就可以忽略， 从而 得到一个近似表达式： 0 )()()()()(rdrurorrgruru bs （2.21） 这个式子就构成了一阶的 born 近似。使用)( )1( rub来表示这个一阶的 born 近似， 同样可以通过将式子（2.21）右边中的)( 0 ru替换成 )1( 0b uu 得到散射场的二阶 born 近似的表达式： )1( 0 )2( )()()()()(rdrururorrgru bb （2.22） 进一步推广到 i+1 阶的情况： )( 0 )1( )()()()()(rdrururorrgru i b i b （2.23） 可以看出，第 i+1 阶的 born 近似由第 i 项构成。高阶 born 近似对散射场的近似 更加精确，但是其在图像重建中的应用研究并不成熟，所以本文中只考虑一阶近 似的情况，简单起见用)( rub表示。 另外，这里特别要注意的是 born 近似只有当散射场的幅度远小于入射场的 情况下才成立。假设物体是一个均匀圆柱体，即其折射系数常数的情况下，可以 将这一条件表示成物体尺寸和折射系数的函数。 假设入射波)( 0 ru是沿方向 0 k传 播的平面波，这样对于一个大型物体，物体内部的场就不能表示为： rkj object aeruru 0 )()( （2.24） 而应该表示成一个反应折射率变化 n的函数。沿着通过圆柱体中心且平行于入 射平面波传播方向的射线，物体内部场就成为入射波变化的情况： rknj object aeru 0 1 )( ）（ （2.25） 由于平面波通过了物体， 因此入射波和物体内部场的相位差就约等于物体内 部折射率变化的积分。对圆柱体来说，总的相移就可以近似为： a nchangphase4 （2.26） 第 2 章 衍射层析成像的基础理论 11 其中 a 是圆柱体半径，是入射波的波数。为了使一阶 born 近似有效的必要条 件就是入射波和穿透物体的波之间的相移小于。这个条件可以表述为： 4 an （2.27） 这个式子反映了入射波波长、折射率变化率以及物体尺寸之间的关系。实际应用 中，我们要选选择合适的波长，使得这三种的关系满足这一关系式。 2.1.3 rytov 近似 rytov 近似成立的条件相对于 born 近似稍微宽松一些。这里我们先将总场 表示为： )( )( r eru （2.28） 并将波动方程重新写成一下形式： nkk 2 0 2 0 22 2)( （2.29） 总相移等于入射波的相移函数 0 同散射场的复相移 s 之和： )()()( 0 rrr s （2.30） 其中 )( 0 0 )( r eru ，然后我们有： 02)(2)( 2 0 2 0 2 0 22 0 2 0 nkk sss （2.31） 和在 born 近似中一样，我们得到以下式子： 0)( 2 00 22 0 k （2.32） 带入式（2.31）得到 nk sss 2 0 22 0 2)(2 （2.33） 这个式子是非齐次且非线性的，但是我们可以通过下面的关系将其变成线性方 程： ssss uuuu 2 000 2 0 2 2)( （2.34） 因为我们有： rkj aeu 0 0 （2.35） 于是得到： s s sss ukuuuu 000 22 000 )(2 （2.36） 将上式带入式（2.33）得到： 第 2 章 衍射层析成像的基础理论 12 2)()( 2 0 2 00 2 0 2 nkuuk ss （2.37） 和之前一样，这个微分方程的解同样可以表示成积分等式： 2 0 2 0 0 2)()( drnkurrgu s v s （2.38） 这里的格林函数由公式（2.14）指定。 根据 rytov 近似，假设上式中括号中的项可以有下列式子近似： nknk s 2 0 2 0 2 22)( （2.39） 这样散射场的相位 s 的一阶 rytov 近似为： 2 00 0 )()( )( 2 )( drnkrurrg ru r v s （2.40） 带入式（2.21） ，得到 ru ru r s s ( )( )( 0 （2.41） 和 born 近似一样，rytov 只有在满足一定的条件下才能成立，我们这里稍 加讨论。相对于 born 近似来说，rotov 近似成立的条件相对宽松一些，在之前 的推导过程中，我们假设： nknk s 2 0 2 0 2 22)( （2.42） 有这个式子只有满足如下关系时成立： 2 0 2 )( k n s （2.43） 取 2 0 k带入则有： 2 2 s n （2.44） 这里和 born 近似不同的是，物体的尺寸不会影响到 rytov。另外，这里是 每个单位波长的相移变化 s 起决定作用， 而不是总的相移。 所以当单位波长上的 相移很小时，rytov 近似是有效。 总的来说，我们看到一阶 born 近似和一阶 rytov 近似的表达式是相似的， 这类表达式也构成了我们重建算法的基础。born 近似中我们测量散射场的复数 幅度并将其作为函数 b u的估计值，而在 rytov 近似中我们将散射场的相位作为 b u的估计值。一般来说，基于相位近似的 rytov 近似的结果会更精确一些。同 第 2 章 衍射层析成像的基础理论 13 时，当小尺寸的物体内部折射系数变换很大时，使用 born 近似将会得到比 ryov 近似更好的结果；相反，当大尺寸物体内部折射系数变换不大时，rytov 近似就 会更加适用；当物体尺寸很大并且其折射系数相对周围介质变换很小时，两种近 似得到的结果是一样的。 2.2 傅里叶衍射傅里叶衍射定理定理 傅里叶散射定理将散射场（投影）的傅里叶变换同物体傅里叶变换域中的圆 弧联系到了一起，构成了我们衍射层析成像的基础。下面我们给出了这个定理的 推导过程。 图图 2.1 平面波入射物体示意图(透射型) 如图 2.1 所示，考虑一个平面波入射到物体，前向散射场可以在接受线处测 得。二维平面波可以表示为： rkj eru)( 0 （2.45） 其中的),( yx kkk 满足关系： 2 0 22 kkk yx （2.46） 我们之前提到二维格林函数的表达式为： )( 4 )|( 00 rrkh j rrg （2.47） 其中 0 h是零阶第一类 hankel 函数。 第 2 章 衍射层析成像的基础理论 14 我们重写散射场的积分等式如下： 0 )()()()(rdrrgrurorub （2.48） 这个等式可以认为是格林函数)( rrg同表示物体的函数)( ro与入射场)( 0 ru 的乘积之间的卷积。我们这里先定义一下傅里叶变换对： )()( )()( )()( kuru kgrrg koro （2.49） 然后波动方程的积分解可以写成傅里叶变换的形式： )()()()( 0 uogus （2.50） 符号代表卷积，其中),( ，波动方程的傅里叶表达式可以表示为： )(2)( 0 ku （2.51） 这样式(2.50)中的卷