不同图像重建算法对定量化光声重建效果的影响

宋宁宁，秦航，蒋红兵，余晶

1.南京医科大学附属南京医院（南京市第一医院） 医疗设备处，江苏 南京 210006；2.南京市卫生信息中心，江苏 南京 210003

引言

光声成像结合了光学成像高对比度和超声成像高分辨度的优点。光声成像的原理是用光照射生物组织，生物组织吸收一部分能量并转化成热能，使生物组织的温度升高，引起热膨胀并产生超声波。生物组织吸收的能量正比于光子数密度和生物组织的吸收系数。光声重建的第一步是用采集到的超声波重建初始场强，初始场强是Grüneisen系数（反映生物组织的声学特性）、光子数密度和生物组织的光吸收系数三者的乘积，只有生物组织的光学系数才能反映生物组织的本质，所以光声成像的第二步是分离出生物组织的光学系数（如光的吸收系数和散射系数）和声学系数（Grüneisen系数等），光声重建的第二步也称定量化光声重建。

对于定量化光声重建，已经有大量的研究成果，我们大体上归为4类：

（1）重建光的吸收系数。假设光的散射系数是已知的，Ripol等[1]用一个点光源重建了光的吸收系数，考虑了光源的长度，光源脉冲的形状，脉冲响应和超声换能器的尺寸有限性。Banerjee等[2]从测量到的边界超声直接重建了光的吸收系数。

（2）重建光的吸收系数和光子数密度。Rosenthal等[3]提取了光的吸收系数和光子数密度，这个算法不是基于光的光扩散方程的具体解，所以对光源的形状没有具体要求。

（3）通过一个不同波长的光源同时重建光的吸收系数和散射系数，Guillaume等[4]已经证明，用一个不同波长的光源可以同时重建光的吸收系数，散射系数和Grüneisen系数。

（4）通过多个光源同时重建光的吸收系数和散射系数。自上世纪以来，学者们做了大量的工作旨在提高光声重建的定量化水平。通过从初始场强的分布中重建出生物组织的不同物理参数：光的吸收系数，光的散射系数和Grüneisen系数。Shao等[5]用多个光源从重建的初始场强中分离出光的吸收系数，光的散射系数和Grüneisen系数；在他们最新的工作中，直接从测量到的超声波中分离出了光的吸收系数和散射系数。

本研究直接从测量到的超声波中分离出了光的吸收系数和散射系数，同时考虑了不同图像重建算法对定量化光声重建效果的影响。

1　光声成像物理理论

下面我们回顾一下光声成像过程中的物理定义和光声成像中的主要公式，包括光声成像的正向过程和逆向过程（光声重建过程）。

在光声成像中，入射光的脉宽必须小于生物组织的热弛豫时间和压力扩散时间，人们把这两个限制条件称之为热限制条件和压力限制条件。只有满足这两个条件，在激光脉冲周期内，生物组织吸收光能产生的热量来不及外散只能用于引起自身的温度升高，同时，温度升高引起的热弹性膨胀压力也来不及外传，只能形成压力波向外传输[6-9]。满足这两个条件的光脉冲时间通常只有几个纳秒，所以，热函数可以写成一个正比于Dirac函数的瞬时函数。在光声成像中，生物组织吸收的热量和光照强度是成比例的。

其中，H(r,t)是热函数，定义为在空间r，时间点t转化成的热能。(r)是光子数密度。μa(r)是吸收系数。光子数密度可通过求解光扩散方程得到。

D(r)是扩散系数，Q(r)是由脉冲激光的照射所产生的光源分布。生物组织吸收的热量转化成热能，组织的温度上升，产生热膨胀，初始场强产生[10-16]：p0(r,u)=Γ(r)μa(r)(r,u)，Γ(r)是Grüneisen系数，表示了生物组织吸收的能量转化为声压的效率[10-16]。

式(3)是外传的超声波的表达式。为了避免计算对超声波对时间的偏导，我们计算的pmodel(r,u,t)是超声波对时间的积分乘以声速的平方。

其中u是生物组织的光学参数。u=[μaD]T。对于有异常的生物组织来说，生物组织的光学参数可以写成u=u0+δu，其中u0是正常组织的光学参数，δu是在正常组织光学参数的基础上叠加的光学参数的异常。如果δu＜＜u0，我们可以把测量到的超声波表达式展开成一阶泰勒级数的形式。

这里Δ=p-pmodel，H是Hessian矩阵。

2　仿真实验及结果

我们以人体乳腺组织为例，进行了下面的仿真实验（图1），选取了图1所示6 cm×6 cm的一个方形测试区域，其中位于中心的2 cm×2 cm的区域用来做光声重建。重建区域正常组织的吸收系数和散射系数分别是μa0=0.1 cm-1和D0=0.03 cm。在正常组织的基础上加入了一个扩散系数的异常δD=0.003 cm（位于中心，尺寸是0.5 cm×0.5 cm），两个吸收系数的异常δμa=0.01 cm-1（尺寸分别是0.3 cm×1.1 cm和0.3 cm×0.3 cm，偏离中心0.6 cm）。选择的异常相对较小（正常组织的10%以内），这样能够满足Born approximation。如果异常相对正常组织较大，需要用非线性重建算法。光源位于重建区域0.3 cm处，超声传感器位于重建区2 cm处，即测试区域的边缘。

我们测试了两种不同的图像重建算法对重建效果的影响：奇异值分解算法（Singular Value Decomposition，SVD）和代数迭代算法（Algebraic Reconstruction Technique，ART）。

在奇异值分解算法中，Hessian矩阵分解成：H=U×S×VT，S是和H同维的对角矩阵，其中非负奇异值iσ是以降序排序的。U和V是单位矩阵。为了避免稀疏矩阵H引起的不稳定性，采用Tikhonov正则化，用代替了S-1的对角值。λ是Tikhonov正则化算子。

代数迭代重建算法用迭代方法寻求δu的解：

式(7)中，λn是一个加速算子，在下面所有的仿真中，λn取常数值1。

图1　测试区域

在下面的仿真实验中，我们用了4个点光源和60个超声传感器，超声传感器均匀的分布在测试区域的边缘（一侧15个传感器）。测试了两种不同的图像重建算法：SVD和ART在不同的噪声水平下对定量化图像重建效果的影响。在不同的噪声水平下重建的吸收系数和散射系数的异常，见图2。

图2　不同算法下重建的光的吸收系数和散射系数

从图2我们可以看出，重建结果对噪声是非常敏感的，不同的噪声水平下重建效果差别很大，在一定噪声水平范围内（本文研究的噪声范围均≤1%），随着噪声的增大，ART的重建效果似乎更好一些，不同噪声水平下重建图像和原图像之间的平方差，见表1。在一定噪声水平范围内，ART的重建效果要比SVD好。在以后的重建中，我们选择ART作为重建算法。

表1　不同算法和噪声水平下重建图像和原图像之间的平方差

3　结论

我们用了多个光源来做定量化光声重建，此方法可以同时重建光的吸收系数和散射系数。在定量化光声重建过程中，我们测试了两种不同的线性图像重建算法：SVD和ART。测试结果表明，将一定噪声水平融进去，ART比SVD效果好。同时，从结果中我们可以看出，这两种线性算法对噪声都是非常敏感的，在噪声进一步加大的情况下，我们应考虑其他非线性算法来提高重建图像质量。

[参考文献]

[1] Ripoll J,Ntziachristos V.Quantitative point source photoacoustic inversion formulas for scattering and absorbing media[J].Phys Rev E,2005,71.

[2] Banerjee B,Bagchi S,Vasu RM.Quantitative photoacoustic tomography from boundary pressure measurements:noniterative recovery of optical absorption coef fi cient from the reconstructed absorbed energy map[J].J Opt Soc Am,2008,25:2347-2356.

[3] Rosenthal A,Razansky D,Ntziachristos V.Quantitative optoacoustic signal extraction using sparse signal represen-tation[J].Med Imaging IEEE Trans,2009,28:1997-2006.

[4] Guillaume B,Ren K.On multi-spectral quantitative photoacoustic tomography in diffusive regime[J].Inverse Probl Imag,2012.

[5] Shao P,Cox B,Zemp RJ.Estimating optical absorption,scattering, and Grueneisen distributions with multipleillumination photoacoustic tomography[J].Appl Opt,2011,50:3145-3154.

[6] 肖嘉莹.定量光声成像技术及在骨关节炎诊断的研究[D].长沙:中南大学,2011.

[7] 向良中,邢达,谷怀民,等.改进的同步迭代算法在光声血管成像中的应用[J].物理学报,2007,56(7):3911-3916.

[8] 吴丹,陶超,刘晓峻.光声成像中延迟求和方法和反重影重构方法的比较[J].无损检测,2011,33(9):36-39.

[9] 杨迪武,黄仲,曾吕明,等.基于环形阵列探测器的快速光声成像[J].激光生物学报,2015,24(5):423-427.

[10] 简晓华,崔崤峣,向永嘉,等.自适应多光谱光声成像技术研究[J].物理学报,2012,61(21):456-461.

[11] 黄盛松,刘博,吴登龙.光声成像在泌尿系疾病的研究进展[J].外科研究与新技术,2016,5(1):49-52.

[12] 关天培,方驰华.光声成像技术及其在原发性肝癌边界界定中的应用[J].中华肝脏外科手术学电子杂志,2016,5(2):65-67.

[13] 吴宁,任秋实,李长辉.光声层析成像研究进展[J].中国医疗设备,2015,3(2):16-20.

[14] 孙正,苑园,王健健.血管内光声成像的研究进展[J].中国生物医学工程学报,2015,(2):221-228.

[15] 翁国星,陈远翔,唐嘉铭,等.光声成像技术评估慢性心肌缺血的活体研究.中国激光医学杂志,2016,25(5):256-259.

[16] Herzog E,Taruttis A,Beziere N,et al.多谱线光声断层摄影术的癌症异质性光学成像[J].中国医疗设备,2012,27(9):21-26.