CN114399562B - 一种结合光通量补偿的迭代量化光声成像方法 - Google Patents

Abstract

本发明公开了一种结合光通量补偿的迭代量化光声成像方法，所述算法包括如下步骤：一、依据光声效应获取光声数据，并进行光声图像重建；二、根据重建图像换算得到光能量分布；三、初始化基于Monte Carlo仿真的模型；四、基于背景的先验知识设置仿真模型参数进行初次光通量补偿；五、更新仿真模型参数；六、基于更新的参数由Monte carlo计算得到光通量；七、计算第k次迭代得到的光能量沉积；八、计算测量得到的光能量沉积与迭代计算得到的光能量之间的误差ε，若误差足够小则停止迭代，输出此时的吸收系数，否则，更新吸收系数，转至步骤五继续迭代。本发明的算法实现了从光能量沉积图中定量恢复出组织中不同吸收体的吸收分布。

Description

一种结合光通量补偿的迭代量化光声成像方法

技术领域

本发明属于生物医学定量光声成像技术领域，涉及一种基于Monte carlo模型的迭代量化组织光学吸收系数分布算法。

背景技术

光声成像是一种基于光学方法激发组织产生超声信号的新型成像方式，具有光学高对比度和超声大穿透深度与高分辨率优势，同时可以获取血红蛋白、脂质成分、血氧代谢等分子和功能信息。此外，光声技术还具有跨尺度成像能力，因而可以灵活实现多种成像模式。光声成像是肿瘤诊断、疗效监测、诊疗一体化等方向从基础到临床研究的关注热点。

量化光声成像技术是目前光声成像领域的研究热点之一。通过光声成像系统直接获得的光声图像本质上是生物组织的光吸收密度图像，它是组织光吸收系数、光通量和Gruneisien参数共同作用的结果。然而，与疾病各种功能和生理参数(如：氧饱和度、血红蛋白浓度、脂肪含量)直接相关的是组织光学系数。此外，借助外源对比剂的分子成像技术是光声成像用于癌症早期诊断的重要工具，通过获取对比剂在组织的浓度分布实现疾病的病灶定位和诊断。因此，通过量化光声成像技术，提取生物组织光学系数图像，是获取生物组织重要功能参数的前提条件，对推进光声成像技术在临床生理功能监测、癌症早期诊断、心脑血管易损斑块检测和分子成像领域的应用具有重要的意义。

发明内容

本发明针对量化过程中光通量非均匀分布问题，结合光通量补偿思想，提供了一种结合光通量补偿的迭代量化光声成像方法。在该迭代算法中，计算每次迭代时被检测到的能量分布与计算得到的吸收能量分布之间的偏差，通过减小吸收能量的偏差，恢复值收敛到真实的吸收分布。

本发明的目的是通过以下技术方案实现的：

一种结合光通量补偿的迭代量化光声成像方法，包括如下步骤：

步骤一、通过阵列光声成像系统依据光声效应获取光声数据，并根据光声图像重建算法进行光声图像重建；

步骤二、根据重建图像换算得到光能量分布H₀，记为Azx0，换算公式为：

PA＝ξ*H；

其中，PA为光声信号，ξ为常数，室温下数值为0.11，H为光能量沉积；

步骤三、根据背景组织的光学特性知识初始化基于Monte Carlo仿真的模型；

步骤四、基于背景的先验知识设置仿真模型参数进行初次光通量补偿，其中：迭代次数设为k次(k≥0)，初始化k＝0，预设吸收系数迭代初始值μa⁰，μa^k＝μa⁰；由Monte carlo计算得到背景组织光通量fluence⁽⁰⁾，记为Fzx⁰，依据

进行初次光通量补偿，计算得到补偿后的吸收系数分布，其中σ为正则化因子；

步骤五、更新仿真模型参数；

步骤六、基于更新的参数由Monte carlo计算得到光通量fluence^(k)，记为Fzx^k；

步骤七、计算第k次迭代得到的光能量沉积H^k，H^k＝μa^k×Fzx^k；

步骤八、计算测量得到的光能量沉积与迭代计算得到的光能量之间的误差ε，若误差足够小(0.005～0.05)则停止迭代，输出此时的吸收系数，否则，更新吸收系数μa^k，转至步骤五继续迭代。

相比于现有技术，本发明具有如下优点：

1、本发明的算法实现了从光能量沉积图中定量恢复出组织中不同吸收体的吸收分布。

2、相比于传统迭代算法，本发明的算法基于Monte Carlo仿真模型保证了定量结果的准确性。

3、本发明的算法改进的参数更新方法使得迭代过程被加速。

附图说明

图1为结合光通量补偿的迭代量化光声成像方法流程图；

图2为光声成像过程示意图；

图3为组织模型z-x切面结构图；

图4为光能量沉积；

图5为背景组织模型；

图6为背景组织对应的光通量分布；

图7为已知待恢复的光能量沉积；

图8为背景组织光通量补偿结果；

图9为吸收系数真实值；

图10为吸收系数恢复值；

图11为恢复吸收系数1、吸收系数2结果；

图12为迭代误差。

具体实施方式

下面结合附图对本发明的技术方案作进一步的说明，但并不局限于此，凡是对本发明技术方案进行修改或者等同替换，而不脱离本发明技术方案的精神和范围，均应涵盖在本发明的保护范围中。

对于采集到的光声信号，可用如下的公式描述：

PA＝Γμ_aΦ (1)；

其中，PA为采集到的光声信号，Γ为格日尼森系数，μ_a为组织的光吸收系数，Φ为对应光通量分布。在一定条件下，格日尼森系数和光通量均为常数，因此光吸收系数与光声信号成比例。结合光声光谱技术可进一步求得对肿瘤研究极为重要的血氧饱和度(SO₂)信息。传统方法在求解血氧饱和度时通常假设光通量为不变量，但实际上光通量随深度变化而变化。本发明中通过对采集到的光声信号进行光通量补偿，矫正光通量变化对精确计算血氧饱和度等功能信息的影响，结合结构量化提供的先验知识，实现更精确的量化肿瘤组织的相关功能参数。

针对非均匀光通量分布对吸收分布重建的限制问题，本发明提出了一种定量重建组织光学吸收系数分布的新方法。该方法基于三维蒙特卡罗模拟，采用迭代算法从光吸收能量图中提取吸收系数，并结合光强补偿得到初始化参数。在迭代算法中，计算每次迭代时被检测到的能量分布与计算得到的吸收能量分布之间的偏差。通过减小吸收能量的偏差，恢复值收敛到真实的吸收分布。从理论上和实验上表明，该方法能够准确定量地估计光吸收系数的分布。

本发明中提到的光能量沉积理论上与光声信号存在比例关系，因此可等同处理，即：

PA＝ξ*H (2)；

其中，PA为光声信号，ξ为常数，室温下数值为0.11，H为光能量沉积。

本发明中，最终要实现的目标是：对已经测得的光声信号(光能量沉积)做光通量补偿，最终恢复出组织的真实的吸收系数。

本发明中，如图1所示，迭代算法的实现步骤如下：

步骤一、通过现有的阵列光声成像系统依据光声效应获取光声数据(由超声换能器采集到的时间序列光声信号)，并根据现有的光声图像重建算法(例如滤波反投影，时间反演算法等)进行光声图像重建。

步骤二、根据重建图像换算得到光能量分布H₀，记为Azx0。

步骤三、根据背景组织的光学特性知识初始化基于Monte Carlo仿真的模型，该模型为基于现有的Monte Carlo仿真平台编程构建组织的仿真模型。

步骤四、基于背景的先验知识设置仿真模型参数进行初次光通量补偿，所述仿真模型参数包括初始化模型参数、迭代次数、吸收系数、散射系数、各向异性因子，其中：迭代次数设为k次(k≥0)，初始化k＝0，预设吸收系数迭代初始值μa⁰(可根据组织中吸收体参数自行修改设置，一般可设置0.5)，μa^k＝μa⁰；由Monte carlo计算得到背景组织光通量fluence⁽⁰⁾，记为Fzx⁰。依据

进行初次光通量补偿，计算得到补偿后的吸收系数分布，其中σ为正则化因子，可设为常数，可在每次迭代之前修改，其作用为使迭代误差收敛。

步骤五、更新仿真模型参数，有如下两种可选方法：

第一种方法是由error＝Azx0-Azx^k，计算测量得到的光能量沉积与本次迭代得到的光能量沉积之间的误差，依据此误差计算迭代步长Δμa，

从而更新吸收系数的值μa^k+1＝μa^k+Δμa。

第二种方法是μa^k+1＝μa^k+Δμa×(k+1)，其中Δμa为迭代步长，k为迭代次数，此为逐步逼近的方法，只要步长足够小，迭代次数足够多，就能准确恢复出吸收系数的值。

此处也可以两种方法相结合使用，可加速迭代过程。通常迭代初期使用第一种方法加速迭代过程，后期使用第二种方法提高准确度并保证收敛，即：

(1)迭代初期，使用下述方法加速迭代过程：

由error＝Azx0-Azx^k，计算测量得到的光能量沉积与本次迭代得到的光能量沉积之间的误差，依据此误差计算迭代步长Δμa，

从而更新吸收系数的值μa^k+1＝μa^k+Δμa；

(2)当误差缩小的趋势逐渐平稳，ΔC1-ΔC2＜0.05时，使用下述方法提高准确度并保证收敛：

μa^k+1＝μa^k+Δμa×(k+1)；其中ΔC1为k-1次迭代误差与k次迭代误差之间的差量，ΔC2为k次迭代误差与k+1次迭代误差之间的差量，Δμa为迭代步长，k为迭代次数。

步骤六、基于更新的参数由Monte carlo计算得到光通量fluence^(k)，记为Fzx^k。

步骤七、计算第k次迭代得到的光能量沉积H^k，H^k＝μa^k×Fzx^k。

步骤八、计算测量得到的光能量分布与迭代计算得到的光能量沉积之间的误差ε，若误差足够小(一般可设置误差阈值为0.005，可根据实际测量目标更改)则停止迭代，输出此时的吸收系数，否则，更新吸收系数μa^k，转至步骤五继续迭代。

本发明中，具体的光声成像理论如下：

光声成像兼具光学成像和超声成像的优点，一方面，生理组织对超声信号的散射要比光信号低2到3个数量级，因此它可以提供较深的成像深度和较高的空间分辨率；另一方面，光声成像根据不同组织对激光的选择性吸收，利用特定波长的激光脉冲对组织进行照射，并间接地对脉冲能量在生理组织中的吸收分布进行成像，因此相比纯超声成像，光声图像中不同组织间的光学对比度较高。光声成像的过程可以大致分为产生光声信号、接收光声信号、处理光声信号三个部分。光声成像原理为光声效应，当短脉冲激光照射生物组织时，生物组织的光吸收体(如黑色素、血红蛋白和水等)吸收一部分激光能量，并将其转换为热能，引起局部温度升高，经过热弹性效应产生以超声波形式发射的光声波。由放置于组织周围的超声换能器接收产生的光声信号，通过信号处理与图像重建形成反映组织内部结构和功能的光声图像。光声成像的具体过程如图2所示。

有效地产生光声信号是重建出质量高的组织光声图像的前提。实现光致超声过程需满足两个限制条件，即热约束和压力约束。组织中的吸收体吸收激光能量转化成热能后，热传导过程中的热扩散时间τ_th可近似为式(3)：

其中，L_p为组织体积的线性特征尺寸；D_T为组织的热扩散系数，对于一般软组织而言D_T＝1.4×10^-3cm²/s。

吸收体吸收一个持续时间为τ_p的激光脉冲后，在该脉冲周期内，热扩散的长度可表示为式(4)：

其中，τ_p为激光脉冲宽度。

根据式(3)和式(4)，热约束可以表述为激光脉冲宽度小于热扩散时间，即τ_p＜τ_th。另外，压力穿过整个受热区域的时间表示为τ_s＝L_p/c，其中c表示声速，则压力约束可以表述为τ_p＜τ_s。

由于热约束和压力约束的存在，一般选择的激光脉宽多为纳米量级。短脉冲的持续时间比吸收体的热扩散和压力扩散的时间要小的多，因此在实际应用中热扩散和压力扩散产生的影响便可被忽略。

当满足以上两个约束条件时，初始声压见式(5)：

P₀(r)＝Γη_thμ_a(r)F (5)；

其中，r为超声换能器到被测位置的距离；Γ为格日尼森系数；η_th为光子能量转化为热量的百分比；μ_a为组织的光吸收系数；F为组织的局部光通量。

初始声压决定了热膨胀过程中某时刻某点的声压值

由式(6)求得：

其中，β为等压体积膨胀系数；

为热函数，表示沉积的热能。

该式忽略了热扩散和运动粘度，等号左边描述了声波的传播过程，右边描述了热能量源。对各个部位各时间的声压信号

进行重建，可获得研究组织的光声图像。

当脉冲激光的脉冲宽度很窄，满足介质中激光能量的沉积时间小于能量扩散时间的是条件，由介质温度变化导致的其他影响可以忽略不计。在理想的脉冲激光的作用下，光声信号的特性由光能量的吸收分布决定，因此初始压力分布P0即光声信号是与吸收能量图即H成比例的。

具体实例：

如下为本发明算法的数值仿真验证，仿真中，可根据H＝μa×Fzx直接计算得到光能量沉积分布。

建模如下：在一个组织块中设置两个大小不同且吸收系数不同的的吸收体，如图3和图4所示。

参数设置如表1所示，两个吸收体具有不同尺寸，不同吸收系数。

表1

tissue list	μacm<sup>-1</sup>	μscm<sup>-1</sup>	g	nnm	sizecm
						small target	2	100	0.8	532	z(0.2--0.3),x(-0.2--0.1)
large target	1	100	0.8	532	z(0.4--0.6),x(0.2--0.4)
						background	0.1	100	0.8	532	1.2*1.2*1.2

注：μa为组织吸收系数，μs为散射系数，g为各向异性因子，n代表波长，size代表尺寸。

测试结果如图5-图12、表2所示。

表2

name	μa1 cm<sup>-1</sup>	μa2 cm<sup>-1</sup>
			迭代恢复值	1.99487	1.02211
真实值	2	1
			恢复误差	1.22％	0.31％
恢复精度	99.74％	97.79％

Claims (4)

1.一种结合光通量补偿的迭代量化光声成像方法，其特征在于所述方法包括如下步骤：

步骤一、通过阵列光声成像系统依据光声效应获取光声数据，并根据光声图像重建算法进行光声图像重建；

步骤二、根据重建图像换算得到光能量沉积H₀，记为Azx0；

步骤三、根据背景组织的光学特性知识初始化基于Monte Carlo仿真的模型；

步骤四、基于背景的先验知识设置仿真模型参数进行初次光通量补偿，其中：迭代次数设为k次，k≥0，初始化k＝0，预设吸收系数迭代初始值μa⁰，μa^k＝μa⁰；由Monte carlo计算得到背景组织光通量fluence⁽⁰⁾，记为Fzx⁰，依据

进行初次光通量补偿，计算得到补偿后的吸收系数分布，其中σ为正则化因子；

步骤五、更新仿真模型参数：

(1)迭代初期，使用下述方法加速迭代过程：

由error＝Azx0-Azx^k，计算测量得到的光能量沉积与本次迭代得到的光能量沉积之间的误差，依据此误差计算迭代步长Δμa，

从而更新吸收系数的值μa^k+1＝μa^k+Δμa；

(2)当误差缩小的趋势逐渐平稳，ΔC1-ΔC2＜0.05时，使用下述方法提高准确度并保证收敛：

μa^k+1＝μa^k+Δμa×(k+1)；其中ΔC1为k-1次迭代误差与k次迭代误差之间的差量，ΔC2为k次迭代误差与k+1次迭代误差之间的差量，Δμa为迭代步长，k为迭代次数；

步骤六、基于更新的参数由Monte carlo计算得到光通量fluence^(k)，记为Fzx^k；

步骤七、计算第k次迭代得到的光能量沉积Azx^k，记为Azx^k＝μa^k×Fzx^k；

步骤八、计算测量得到的光能量沉积与迭代计算得到的光能量沉积之间的误差error，若误差小于或等于误差阈值则停止迭代，输出此时的吸收系数，否则，更新吸收系数μa^k，转至步骤五继续迭代。

2.根据权利要求1所述的结合光通量补偿的迭代量化光声成像方法，其特征在于所述步骤一中，光声图像重建算法为滤波反投影算法或时间反演算法。

3.根据权利要求1所述的结合光通量补偿的迭代量化光声成像方法，其特征在于所述步骤二中，换算公式为：

PA＝ξ*H；

其中，PA为光声信号，ξ为常数，H为光能量沉积。

4.根据权利要求1所述的结合光通量补偿的迭代量化光声成像方法，其特征在于所述误差阈值为0.005。

Images