一种用于生物光声内窥成像的光吸收系数重建方法

摘要

一种用于生物光声内窥成像的光吸收系数重建方法，属于医学成像技术领域，其技术方案是，所述方法首先根据超声探测器从各个角度采集的生物腔体组织产生的光声信号重建出光吸收能量分布的测量值；然后通过前向仿真获得光吸收能量分布的理论值；最后根据光吸收能量分布的测量值和理论值构造目标函数，并通过对目标函数进行优化，重建出腔体横截面上组织光吸收系数的空间分布。本发明能够利用生物腔体组织产生的光声信号完整地重建出腔体横截面上组织光吸收系数的空间分布，为疾病的早期诊断提供机体组织变化的准确可靠的参考信息。

权利要求

1.一种用于生物光声内窥成像的光吸收系数重建方法，其特征是，所述方法首先根据采集的光声信号重建出光吸收能量分布的测量值；然后通过前向仿真获得光吸收能量分布的理论值；最后根据光吸收能量分布的测量值和理论值构造目标函数，并通过对目标函数进行优化，重建出腔体横截面上组织光吸收系数的空间分布；

所述方法按以下步骤进行：

a.重建光吸收能量分布的测量值

在直线扫描模式下，根据超声探测器从各个测量角度采集到的光声信号，得到光吸收能量分布的测量值：其中，r是θ-l平面极坐标系中的一点，c是光在组织内的传播速度，Hm(r)是位置r处的光吸收能量分布的测量值，z0是位置r处与组织表面θ轴之间的垂直距离，ri是与角度θi(i＝1,

2,...,m)相对应的位置，cs是超声信号在组织中传播速度，β是组织的等压膨胀系数，Cp是组织的比热容，pi(r,t)是超声探测器在时刻t、角度θi(i＝1,2,...,m)、位置r处采集到的组织产生的光声信号的声压；

b.通过前向仿真获得光吸收能量分布的理论值

首先，用由节点表示的面积元对成像区域进行离散化，并根据成像区域内每个节点处光扩散方程的差分形式，计算出位置r处的光子密度函数Φ(r)；

然后根据位置r处的光子密度函数Φ(r)得到光吸收能量密度的理论值H(r,μa(r))：H(r,μa(r))＝μa(r)·h·f·Φ(r)

其中，μa(r)是位置r处组织的光吸收系数，h是普朗克常数，f是入射光的频率；

c.重建光吸收系数的空间分布

假定组织的光散射系数已知，利用下式求得光吸收系数的空间分布：其中，C是与光声信号采集系统有关的标定因子， 是估算出的、待测区域内位置r处的光吸收系数。

说明书

一种用于生物光声内窥成像的光吸收系数重建方法

技术领域

[0001] 本发明涉及一种利用生物腔体组织产生的光声信号重建组织光吸收系数空间分布的方法，属于医学成像技术领域。

背景技术

[0002] 光声层析成像是一种新型的非电离式生物医学功能成像方法，它结合了超声成像的高分辨率和光学成像的高对比度的优势。PAT成像方法的逆问题包括声学和光学两方面：声学逆问题是指根据超声探测器采集到的光声信号(本质是超声波)重建组织内部的初始声压分布或空间光吸收能量密度；光学逆问题是指根据超声探测器采集到的光声信号或者据此重建出的光吸收能量密度，估算组织光学特性参数的空间分布，得到对组织光学特性的定量评价。光吸收能量密度是由局部的光吸收系数和光子数分布共同决定的，并不能反映组织的光学特性，而光吸收系数与组织的化学成分密切相关，正常和病变组织的光学特性参数之间通常有较明显的差异，因此组织光吸收系数的空间分布图可为疾病的早期诊断提供更加准确可靠的基础参考信息。但到目前为止，组织光吸收系数的重建方法仍不成熟，还有必要进一步进行研究。

发明内容

[0003] 本发明的目的在于针对现有技术之弊端，提供一种用于生物光声内窥成像的光吸收系数重建方法，以获取不同组织的光吸收系数之间的差异，为疾病的早期诊断提供机体组织变化的准确可靠的参考信息。

[0004] 本发明所述问题是以下述技术方案解决的：

[0005] 一种用于生物光声内窥成像的光吸收系数重建方法，所述方法首先根据超声探测器从各个角度采集的生物腔体组织产生的光声信号重建出光吸收能量分布的测量值；然后通过前向仿真获得光吸收能量分布的理论值；最后根据光吸收能量分布的测量值和理论值构造目标函数，并通过对目标函数进行优化，重建出腔体横截面上组织光吸收系数的空间分布。

[0006] 上述用于生物光声内窥成像的光吸收系数重建方法，所述重建按以下步骤进行：

[0007] a.重建光吸收能量分布的测量值

[0008] 在直线扫描模式下，根据超声探测器从各个测量角度采集到的周围组织产生的光声信号，得到光吸收能量分布的测量值：

[0009]

[0010] 其中，r是θ-l平面极坐标系中的一点，c是光在组织内的传播速度，Hm(r)是位置r处的光吸收能量分布的测量值，z0是位置r处与组织表面θ轴之间的垂直距离，ri是与角度θi(i＝1,2,...,m)相对应的位置，cs是超声信号在组织中传播速度，β是组织的等压膨胀系数，Cp是组织的比热容，pi(r,t)是超声探测器在时刻t、角度θi(i＝1,2,...,m)、位置r处采集到的组织产生的光声信号的声压；

[0011] b.通过前向仿真获得光吸收能量分布的理论值

[0012] 首先，用由节点表示的面积元对成像区域进行离散化，并根据成像区域内每个节点处光扩散方程的差分形式，计算出位置r处的光子密度函数Φ(r)；

[0013] 然后根据位置r处的光子密度函数Φ(r)得到光吸收能量密度的理论值H(r,μa(r))：

[0014] H(r,μa(r))＝μa(r)·h·f·Φ(r)

[0015] 其中，μa(r)是位置r处组织的光吸收系数，h是普朗克常数，f是入射光的频率；

[0016] c.重建光吸收系数的空间分布

[0017] 假定组织的光散射系数已知，利用下式求得光吸收系数的空间分布：

[0018]

[0019] 其中，C是与光声信号采集系统有关的标定因子， 是估算出的、待测区域内位置r处的光吸收系数。

[0020] 本发明能够利用生物腔体组织产生的光声信号完整地重建出腔体横截面上组织光吸收系数的空间分布，为疾病的早期诊断提供机体组织变化的准确可靠的参考信息。

附图说明

[0021] 下面结合附图对本发明作进一步说明。

[0022] 图1是含有脂质斑块的血管横截面示意图；

[0023] 图2是θ-l极坐标平面内网格划分的内点示意图；

[0024] 图3是θ-l极坐标平面内网格划分的边界点示意图；

[0025] 图4是腔体成像区域内面向光源的边界点示意图。

[0026] 文中各符号为：θ、极角；l、极径；θ-l、平面极坐标系；m、血管横截面被等角度分割的总份数；θi、成像导管的第i个测量角度，其中i＝1,2,...,m；r、θ-l平面极坐标系中的一点；Hm(r)、位置r处的光吸收能量分布的测量值；z0、位置r处与组织表面θ轴之间的垂直距离；ri、θ-l平面中与角度θi(i＝1,2,...,m)相对应的位置；cs、超声信号在生物软组织中的传播速度；β、组织的等压膨胀系数；Cp、组织的比热容；pi(r,t)、超声探测器在时刻t、角度θi、位置r处采集到的组织产生的光声信号的声压，其中i＝1,2,...,m； 哈密顿算子；Φ(r)、位置r处的光子密度函数；Φ(r,rs,q)、时刻t位置r处的时域光子密度函数的Laplace变换；δ(r-rs)、位置r处无向点源的Laplace变换；q、Laplace变换的复频率因子，当q＝0时，DE是稳态扩散方程；Ω、待成像的组织区域；c、光在组织内的传播速度； 组织表面 点处的外法向量；Qs、光源强度；rs、θ-l平面极坐标系中光源所在位置；Rf、扩散传输内反射系数；n、组织相对于环境的相对折射率；μa(r)、位置r处组织的光吸收系数；μs(r)、位置r处组织的光散射系数；μ′s(r)、位置r处组织的约化散射系数；g、组织的各向异性因子；D、Ω中一个充分光滑的区域；hθ、沿θ轴的正向步长；hl、沿l轴的正向步长；P、区域D中的一点；(ihθ,jhl)、网格划分后θ-l坐标系中点P的坐标；Φi,j、节点(ihθ,jhl)处的光子密度函数值；ki+1,j、点((i+1)hθ,jhl)处的k值；ki,j+1、点(ihθ,(j+1)hl)处的k值；ki-1,j、点((i-1)hθ,jhl)处的k值；ki,j-1、点(ihθ,(j-1)hl)处的k值；|AB|、线段AB的长度；|BE|、线段BE的长度；|EA|、线段EA的长度；|D|、区域D的面积；H(r,μa(r))、位置r处光吸收能量密度的理论值；h、普朗克常数；f、入射光的频率；C、与光声信号采集系统有关的标定因子； 待测区域内位置r处估算出的光吸收系数；X、待求的μa(r)；f(X)、待优化的目标函数； 位置r处光吸收系数的估计值；λ、正则化参数；X0、p0、光吸收系数和次梯度的初始值；Xn+1、pn+1、n次迭代后的光吸收系数和次梯度；Xn、pn、n-1次迭代后的光吸收系数和次梯度；、pn与X的内积；M、Ne×N维的稀疏矩阵；Ne、迭代次数；N、网格个数；|·|、L1准则；α、惩罚参数；ν、第n次迭代产生的残余量，其初始值v0＝0；shrink(·)、shrink算子；wn、第n次搜索的方向；B0、近似逆Hessian矩阵的初始值；I、单位矩阵；Bn、第n次迭代的近似逆Hessian矩阵；an、步长； g(X)在Xn的梯度。

具体实施方式

[0027] 本发明的处理步骤包括：

[0028] 1.由超声探测器接收到的光声信号重建光吸收能量分布的测量值

[0029] 如附图1所示，以血管内光声成像为例，成像导管(超声探测器位于成像导管顶端，用于接收周围组织产生的光声信号)位于血管横截面的中心，周围由内向外依次为管腔、粥样硬化斑块、血管壁内膜/中膜和外膜。忽略超声探测器的孔径效应，将其看作理想的点换能器，其扫描轨迹为平行于成像平面的圆形轨迹。

[0030] 血管横截面所在的坐标系是θ-l极坐标系，其中θ是极角，l是极径，坐标原点是成像导管中心，水平向右的方向为θ轴正方向。以坐标原点为中心，将血管横截面等角度分成m个扇区，成像导管的第i个测量角度是θi＝360(i-1)/m，其中i＝1,2,...,m。

[0031] 在直线扫描模式下，根据超声探测器采集到的光声信号，得到光吸收能量分布的测量值：

[0032]

[0033] 其中，r是θ-l平面极坐标系中的一点，c是光在组织内的传播速度，Hm(r)是位置r处的光吸收能量分布的测量值，z0是位置r处与组织表面θ轴之间的垂直距离，ri是与角度θi(i＝1,2,...,m)相对应的位置，cs是超声信号在生物软组织中的传播速度，β是组织的等压膨胀系数，Cp是组织的比热容，pi(r,t)是超声探测器在时刻t、角度θi(i＝1,2,...,m)、位置r处采集到的组织产生的光声信号的声压。

[0034] 2.通过前向仿真获得光吸收能量分布的理论值

[0035] 首先，用由节点表示的面积元对成像区域进行离散化，并根据成像区域内每个节点处光扩散方程的差分形式，计算出位置r处的光子密度函数Φ(r)。具体步骤如下：

[0036] 采用准直光源模型，将光源视为组织边界处的内向光子流，进而嵌入到Robin边界条件中，则边界含有光源项的DE复频域表达式为：

[0037]

[0038] 其中， 是哈密顿算子；r是θ-l坐标系中的一点；Φ(r,rs,q)是时刻t位置r处的时域光子密度函数的Laplace变换；δ(r-rs)是位置r处无向点源的Laplace变换；q是Laplace变换的复频率因子，当q＝0时，DE是稳态扩散方程；Ω是待成像的组织区域；c是光在组织内的传播速度；是组织表面 点处的外法向量；Qs是光源强度；rs是θ-l平面极坐标系中光源所在位置；Rf是扩散传输内反射系数：

[0039] Rf≈-1.4399n-2+0.7099n-1+0.6681+0.0636n  (3)

[0040] 式(3)中，n是组织相对于环境的相对折射率；μa(r)和μs(r)分别是位置r处组织的光吸收系数和散射系数；当μs'＞＞μa时

[0041]

[0042] 其中，μs′(r)是位置r处组织的约化散射系数：

[0043] μ′s(r)＝μs(r)(1-g)  (5)

[0044] 其中g是组织的各向异性因子。

[0045] 在Ω中取一个充分光滑的区域D∈Ω，在θ-l平面内，对区域D进行矩形网格划分，沿θ轴和l轴的正向步长分别为hθ和hl，D中的一点P在θ-l坐标系中的坐标是(ihθ,jhl)，在区域D上，对式(2)进行积分，得到：

[0046]

[0047] 如附图2所示，阴影部分为区域D，若P是内点，则式(6)的差分形式为：

[0048]

[0049] 其中，Φi,j是节点(ihθ,jhl)处的光子密度函数值；根据式(4)求出ki,j，其中，ki+1,j是在点((i+1)hθ,jhl)处的k值，ki,j+1是在点(ihθ,(j+1)hl)处的k值，ki-1,j是在点((i-1)hθ,jhl)处的k值，ki,j-1是在点(ihθ,(j-1)hl)处的k值。

[0050] 如附图3所示，由弧线 线段AB和BE组成区域D，若P是边界点，则式(6)的差分形式为：

[0051]

[0052] 式中，|AB|、|BE|和|EA|分别是线段AB、BE和EA的长度；|D|是区域D的面积。

[0053] 如附图4所示，光源沿l轴即径向入射，由弧线 线段AB和BE组成区域D，若P是面向光源的边界点，则式(6)的差分形式为：

[0054]

[0055] 然后，根据位置r处的光子密度函数Φ(r)得到光吸收能量密度的理论值H(r,μa(r))：

[0056] H(r,μa(r))＝μa(r)·h·f·Φ(r)  (10)

[0057] 其中，h是普朗克常数，f是入射光的频率。

[0058] 3.重建光吸收系数的空间分布

[0059] 假定组织的光散射系数已知，求解光吸收系数分布的问题表述为如下的非线性最小二乘问题：

[0060]

[0061] 其中，C是与光声信号采集系统有关的标定因子， 是估算出的、待测区域内位置r处的光吸收系数。为表述方便，用X表示μa(r)，并定义：

[0062] f(X)＝||Hm(r)-C·H(r,X)||2  (12)

[0063] 该最优化问题是一个病态问题，本发明方法采用TV正则化(Gao  H,Zhao H.Multilevel bioluminescence tomography based on radiative transfer equation Part 1:l1 regularization.Optics Express,2010,18(3):1854-1871.)使其良态化，将待优化的目标函数改写为：

[0064]

[0065] 其中， 是位置r处光吸收系数的估计值；λ是正则化参数，在迭代过程中设定λ＝1。

[0066] 本发明采用Bregman方法(Gao H,Zhao H,Osher S.Bregman methods in quantitative photoacoustic tomography.Cam Report,2010,30(6):3043-3054)迭代求解式(13)的最优解。具体步骤如下：

[0067] 设定各参数的初始值：光吸收系数X0＝0，次梯度p0＝0。第n次迭代后，光吸收系数Xn+1与次梯度pn+1的结果为：

[0068]

[0069]

[0070] 其中，是次梯度pn与X的内积。

[0071] 在离散网格条件下，式(14)右端的TV项简化为

[0072] ||X||TV＝|MX|  (16)

[0073] 其中，M是Ne×N维的稀疏矩阵，Ne是迭代次数，N是网格个数，|·|是L1准则。

[0074] 令

[0075] dn＝MXn  (17)

[0076] 将式(14)转化为无约束问题：

[0077]

[0078] 其中，α是惩罚参数，取值为1。对式(18)进行n次迭代后的结果为：

[0079]

[0080]

[0081] νn+1＝νn+dn+1-MXn+1(21)

[0082] 其中，vn+1是第n次迭代产生的残余量，其初始值ν0＝0；shrink算子的定义为：

[0083]

[0084] 设

[0085] 求解Xn+1＝argmin[g(X)]最小化问题的具体步骤如下：

[0086] 步骤1：初始化，设X0＝0，近似逆Hessian矩阵的初始值B0＝I；

[0087] 步骤2：计算搜索方向： 其中Bn是第n次迭代的近似逆Hessian矩阵；

[0088] 步骤3：沿搜索方向找到下一个迭代点：

[0089]

[0090] 其中，an是步长； 是g(X)在Xn的梯度：

[0091]

[0092] 的更新方法是：

[0093]

[0094] 步骤4：根据 判断迭代是否停止；

[0095] 步骤5：更新搜索方向wn。