CN101214156B - 声速不均匀介质热声成像的重建算法 - Google Patents

Abstract

本发明属于热声成像技术领域，具体为一种声速不均匀介质热声成像的图像重建方法。该方法通过对声速均匀介质热声成像的时域重建法进行修正，使声波在介质中两点间的传播时间不再与距离成正比，而是通过估计出的声速分布计算出来。当前对声速不均匀介质热声成像的重建算法大都需要声速分布的先验知识，并采用迭代方法计算，速度很慢。本发明无需介质内声速分布的先验知识，而且是一种直接的反向方法，可以一次性快速重建出声速不均匀介质的整幅热声图像。

Description

声速不均匀介质热声成像的重建算法

技术领域

本发明属于热声成像技术领域，具体为一种声速不均匀介质热声成像的图像重建算法。

技术背景

热声成像(往往也称为光声成像)结合了光学成像高对比度和超声成像高分辨率的优点，近年来发展迅速，成为一种广泛关注的新型无损医学成像技术。热声成像非常适用于肿瘤检测、血管成像等方向，也被尝试应用于老鼠脑部的结构和功能成像、分子成像、流速检测等方向。

热声成像中，生物组织受电磁波(常用微波或激光)脉冲照射后产生热膨胀，从而激发出热声波。热声波是一种超声波，由超声换能器探测后，可用来重建生物组织内部的电磁波吸收分布图像。热声成像技术的关键是成像方法，即如何由探测到的热声波重建出组织的电磁波吸收分布。目前已提出很多成像的重建算法，但是这些算法几乎都假定生物组织内部的声速是均匀的。事实上，这个假设有时是不成立的，例如对乳房成像时声速差异可达10％。此时若假定声速均匀来进行图像重建，将会导致图像内目标的错位和模糊^[1，2]。

近两年已提出补偿声速不均匀性的热声成像算法^[1-4]。这些算法大都需要声速分布的先验知识^[1-3]，而实际应用中往往是不知道声速分布的。目前仅有的无需声速分布先验知识的方法^[4]是基于有限元法、通过迭代方法进行计算，重建速度很慢。本发明不仅无需介质内声速分布的先验知识，而且是一种直接的反向方法，可以一次性快速重建出声速不均匀介质的整幅热声图像。

发明内容

本发明的目的在于提出一种无需介质内声速颁布先验知识，速度快的声速不均匀介质热声成像的重建算法。

本发明提出的一种声速不均匀介质热声成像的重建方法。其具体步骤为：对声速均匀介质热声成像的时域重建法进行修正，使声波在介质中两点间的传播时间不再与距离成正比，而是一个待定值；估计出介质中的声速分布，通过微扰理论计算出声波在介质中两点间的传播时间，然后进行图像重建。

下面对各步骤作进一步具体描述。

通常电磁波脉冲持续时间远短于组织的热扩散时间，因此可忽略热扩散。下面给出三维下热声成像的基本方程^[2，5]：

${\▿}^{2}p(r,t)-\frac{1}{c{\left(r\right)}^2}\frac{{\∂}^{2}p(r,t)}{{\∂t}^2}=-\frac{\mathrm{\β}}{C_p}A\left(r\right)\·\frac{\∂I\left(t\right)}{\∂t}---\left(1\right)$

其中p(r，t)是位置r处的声压，A(r)是电磁波吸收分布，t是时间，I(t)是电磁波能量脉冲函数，β是等压膨胀系数，c(r)是声速分布，C_p是样品组织的比热。热声成像算法是一个典型的逆问题，即：如何由p(r，t)求出A(r)。

当超声换能器沿圆球面轨迹r₀(球心为原点，半径为r₀)扫描时，p(r，t)中r＝r₀。通常假定电磁波脉冲函数是一个δ函数，即I(t)＝δ(t)；声速是一个常数，即c(r)＝c。此时，热声成像的时域重建法为^[5]：

$A\left(r\right)=-\mathrm{\η}\underset{{\mathrm{\Ω}}_0}{\∫}\∫d{\mathrm{\Ω}}_0\frac{1}{t}\frac{\∂p(r_0,t)}{\∂t}{|}_{t=|r_0-r|/c_0}---\left(2\right)$

其中η是一个常数，Ω₀是积分面。

事实上，(2)式的物理意义非常明确。位置r处产生的热声波到达超声换能器位置r₀需要的时间为t。因此A(r)等于所有探测位置的热声波在对应时间t上的值(1/t)(

p(r₀，t)/

t)进行投影。如果组织中声速是均匀的，t和r与r₀间的距离成正比。如果声速不是均匀的，t是关于r与r₀的函数，即t＝T(r，r₀)。因此(2)式被修正为：

$A\left(r\right)\≈-\mathrm{\η}\underset{{\mathrm{\Ω}}_0}{\∫}\∫d{\mathrm{\Ω}}_0\frac{1}{t}\frac{\∂p(r_0,t)}{\∂t}{|}_{t=T(r,r_0)}---\left(3\right)$

(3)式就是本发明的反向重建公式，它考虑到了热声成像中的声速不均匀性。那么现在的问题是如何求出T(r，r₀)。

生物软组织内的声速差异一般不大，因此组织的慢度u(r)＝l/c(r)可以写成：

u(r)＝u₀+ε·u₁(r)    (4)

其中u₀是常数，ε是一个很小的值。根据微扰理论，T(r，r₀)可以写成微扰级数的形式^[6，7]：

T(r，r₀)＝T₀(r，r₀)+ε·T₁(r，r₀)+ε²·T₂(r，r₀)    (5)

此处忽略了二阶以上的小量。可以进一步求出^[7]：

$T_0=\underset{l}{\∫}{u}_0\mathrm{ds},T_1=\underset{l}{\∫}{u}_1\mathrm{ds},T_2=\underset{l}{\∫}\frac{1}{2u_0}(u_1^2-{|\▿T_1|}^2)\mathrm{ds}---\left(6\right)$

其中l是r和r₀间的线段，即声速均匀时的声波传播路径。

从(6)式可见，T₀是没有慢度扰动u₁时的声波传播时间；T₁是一阶时间扰动，表示了u₁沿l的影响；T₂是二阶时间扰动，表示了声波折射的主要影响。显然T₂的计算较复杂，因此在生物软组织的热声成像中，考虑到声速差异较小，可以忽略T₂的影响，通过T₀和T₁计算T就足够了。但当介质中的声速差异较大时，例如对生物体的脑部成像时，就需要考虑T₂的影响。

(6)式需要知道慢度(或声速)分布U(r)(或c(r))，而实际应用中往往不知道u(r)。因此本发明提出了一种估计u(r)的方法，从而进一步计算出T(r，r₀)。具体步骤如下。

超声换能器接收到的声波可写为^[2]：

$p(r_0,t)={\mathrm{\η}}_2\frac{\∂}{\∂t}\underset{{\mathrm{\Ω}}_0}{\∫}\∫\frac{A\left(r\right)}{t}d{\mathrm{\Ω}}_0{|}_{T(r,r_0)=t}---\left(7\right)$

其中η₂是一个常数。定义：

$S(r_0,t)=\left[{\∫}_0^{t}p(r_0,t)\mathrm{dt}\right]\·t={\mathrm{\η}}_2\underset{{\mathrm{\Ω}}_0}{\∫}\∫A\left(r\right)d{\mathrm{\Ω}}_0{|}_{T(r,r_0)=t}---\left(8\right)$

可见S(r₀，t)是A(r)的面积分，积分面Ω₀上所有点到r₀的声波传播时间为t。若c(r)是常数，则积分面Ω₀。可简化为球面。如图1所示，当待测组织远比r和r₀间距离小时，Ω₀可近似为一个平面。那么显然有：

S(r₀，t)≈S(-r₀，t′)              (9)

其中t+t′-T(-r₀，r₀)。若声速有较小差异，或者待测组织较大时，(9)仍近似成立。则T(-r₀，r₀)可通过下式计算：

$T(-r_0,r_0)=\arg \max \left(R_{r_0}\left(t\right)\right)---\left(10\right)$

其中R_r0(t)是S(r₀，t)和S(-r₀，-t)的相关函数，定义为：

$R_{r_0}\left(t\right)={\∫}_{-\∞}^{\∞}S(r_0,\mathrm{\τ})S(-r_0,-(\mathrm{\τ}-t))\mathrm{d\τ}---\left(11\right)$

可以认为-r₀与r₀间的声速是近似均匀的，则有：

u_e(r)＝T(-r/|r|·r₀，r/|r|·r₀)/2r₀          (12)

其中u_e(r)表示估计的慢度分布u(r)。

综上，结合(3)、(6)、(10)、(12)式，即可完成声速不均匀的图像重建。首先根据(10)、(12)式用探测到的热声波估计出慢度(或声速)分布；然后根据(6)式计算组织内两点间的声波传播时间；最后根据(3)式重建出待测组织的图像。实际应用中参数ε可随意选取，其结果都一样；参数η、η₂应根据待测组织的等压膨胀系数、比热、声速以及扫描半径来确定^[5]。

附图说明

图1、r和r₀位置处接收到的热声波所对应的待测组织部分。

图2、待测组织的电磁波吸收分布(a)和声速分布(b)。

图3、使用不同声速的图像重建结果，(a)c＝1500m/s，(b)c＝1429m/s，(c)c＝1364m/s，(d)已知声速分布。

图4、声速分布未知时，估计的声速分布(a)和重建的电磁波吸收分布(b)。

具体实施方式

为验证本发明的方法，在计算机上进行二维情况下的仿真实验。虽然本发明的理论论证是在三维情况下进行的，但将其应用到二维下显然也是成立的。仿真实验的具体步骤和每步的结果如下所述：

1、建立待测样品的模型。图2显示了待测组织的电磁波吸收分布和声速分布(m/s)。该声速分布类似于女性乳房模型^[1]。图2的尺寸为18mm×18mm。待测组织的其余参数(如等压膨胀系数、比热)对于成像算法的影响体现在常数η、η₂上，其选取可以不用与真实情况一致，此处简单设置为1。

2、利用时域有限差分法仿真出扫描圆周上接收到的声波，扫描圆周半径为9mm，共有160个等间距的探测位置。

3、先假设声速分布是已知的，用仿真的声波重建出待测样品的图像。图3显示了用均匀声速和已知不均匀声速分布的图像重建结果。

4、声速分布未知时，进行图像重建。先利用本发明的(10)式和(12)式来估计声速分布，然后再基于(3)式和(6)式进行图像重建。估计的声速分布和重建的电磁波吸收分布图像如图4所示。

由步骤3的结果可见，图3(a)～(c)使用均匀声速，结果不好。其中以图3(b)效果相对最好，但最小的圆形组织已经模糊。图3(d)使用本发明的(3)式和(6)式进行图像重建，效果很好。这说明本发明的修正的时域重建法是有效的。重建时，仅对扫描圆周内的区域进行图像重建，而扫描圆周外在重建时像素直接设为0。由于扫描圆周内区域在重建时存在误差，像素本应为0的背景区域重建出的结果是很小的负值。因此图3边缘处有一圆形轮廓，这仅表明扫描轨迹，并不是反映待测物体的轮廓。

由步骤4的结果可见，图4(a)中估计的声速分布与实际情况差别较大，但用在重建电磁波吸收分布图像时，可以较准确地计算出(6)式，获得满意的图像重建结果。这与本发明的理论相符。图4(b)中组织轮廓较清晰，最小的圆形组织也基本没有模糊。总的来看，图4(b)的结果略差于图3(d)，好于图3(a)～(c)。这是因为图4(b)是基于估计出的声速分布进行重建，肯定不及已知声速分布的图3(d)，但会优于对声速不均匀性不进行任何补偿的图3(a)～(c)。

根据实验结果可见，本发明的方法无需声速分布的先验知识，而且可以一次性快速重建出声速不均匀介质的整幅热声图像，效果很好。本发明的仿真实验中声速差异约为10％，一般生物软组织的声速差异会更小一些，因此本方法非常适用于生物软组织的热声成像。

参考文献

[1]Y.Xu，and L.V.Wang，“Effects of acoustic heterogeneity in breast thermoacoustic tomo-graphy，”IEEE Trans.Ultrason.，Ferroelect.，Freq.Contr.，vol.50，no.9，Sep.2003.

[2]X.Jin，and L.V.Wang，“Thermoacoustic tomography with correction for acoustic speedvariations，”Phys.Med.Biol.，vol.51，pp.6437-6448，2006.

[3]J.Zhang，and M.A.Anastasio，“Reconstruction of speed-of-sound and electromagnetic ab-sorption distributions in photoacoustic tomography，”in Proc.SPIE，vol.6080，pp.608619-1-7，2006.

[4]H.Jiang，Z.Yuan，and X.Gu，“Spatially varying optical and acoustic property reconstructionusing finite-element-based photoacoustic tomography，”J.Opt.Soc.Am.，vol.23，no.4，pp.878-888，Apr.2006.

[5]M.Xu，Y.Xu，and L.V. Wang，“Time-domain reconstruction algorithms and numericalsimulations for thermoacoustic tomography in various geometries，”IEEE Trans.Biomed.Eng.，vol.50，no.9，pp.1086-1099，Sep.2003.

[6]R.Snieder，and M.Sambridge，“Ray perturbation theory for traveltimes and ray paths in 3-Dheterogeneous media，”Geophys.J.Int.，vol.109，pp.294-322，1992.

[7]R.Snieder，and D.F.Aldridge，“Perturbation theory for travel times，”J.Acoust.Soc.Am.，vol.98，no.3，pp.1565-1569，Sep.1995.

Claims (1)

1.一种声速不均匀介质热声成像的图像重建方法，其特征在于：对声速均匀介质热声成像的时域重建法进行修正，使声波在介质中两点间的传播时间不再与距离成正比，而是一个待定值；估计出待测介质中的声速分布，通过微扰理论计算出声波在介质中两点间的传播时间，然后进行图像重建；其中：

所述的修正的时域重建法为：

$A\left(r\right)\≈-\mathrm{\η}\underset{{\mathrm{\Ω}}_0}{\∫}\∫d{\mathrm{\Ω}}_0\frac{1}{t}\frac{\∂p(r_0,t)}{\∂t}{|}_{t=T\left({r,r}_0\right)}---\left(3\right)$

这里，A(r)是介质的电磁波吸收分布，r₀是以原点为圆心、半径为r₀的圆球面轨迹，p(r₀，t)是在圆球面上探测到的热声波，t是时间，T(r，r₀)是两点r与r₀间的声波传播时间，η是一个常数；

当待测介质的声速c(r)不均匀时，有：

u(r)＝u₀+ε·u₁(r)                                   (4)

这里，u(r)＝1/c(r)为介质的慢度，u₀是常数，ε是一个很小的值；由微扰理论计算T(r，r₀)的方法为：

T(r，r₀)＝T₀(r，r₀)+ε·T₁(r，r₀)+ε²·T₂(r，r₀)            (5)

$T_0=\underset{l}{\∫}{u}_0\mathrm{ds},$ $T_1=\underset{l}{\∫}{u}_1\mathrm{ds},$ $T_2=\underset{l}{\∫}\frac{1}{2u_0}(u_1^2-{|\▿T_1|}^2)\mathrm{ds}---\left(6\right)$

其中l是r和r₀间的线段；

估计介质的慢度分布u_e(r)由下式估计：

u_e(r)＝T(-r/|r|·r₀，r/|r|·r₀)/2r₀                   (12)

其中计算T(-r₀，r₀)的方法为：

$T(-r_0,r_0)=\arg \max \left(R_{r_0}\left(t\right)\right)---\left(10\right)$

$R_{r_0}\left(t\right)={\∫}_{-\∞}^{\∞}S(r_0,\mathrm{\τ})S(-r_0,-(\mathrm{\τ}-t))\mathrm{d\τ}---\left(11\right)$

是S(r₀，t)和S(-r₀，-t)的相关函数，S(r₀，t)定义为：

$S(r_0,t)=\left[{\∫}_0^{t}p(r_0,t)\mathrm{dt}\right]\·t---\left(8\right).$

Images