一种血管内光声图像的时间反演重建方法与流程

本发明涉及一种对血管内光声图像进行时间反演重建，得到血管壁的轴向截面光吸收分布灰阶图像的方法，属于医学成像技术领域。

背景技术：

血管内光声(intravascular photoacoustic，IVPA)成像技术是继血管内超声(intravascular ultrasound，IVUS)和血管内光学相干断层成像(intravascular optical coherence tomography，IV-OCT)之后的一种新兴微创介入血管成像方法，它兼具纯光学成像和纯超声成像的优点，弥补了现有介入血管成像技术的不足，可以实现血管腔、管壁和易损斑块的高分辨力和高对比度的深度成像。

IVPA的成像原理是将微型光声成像导管插入待测血管腔内，当导管绕其轴旋转时，脉冲激光或幅度调制的连续激光均匀照射到血管壁上，激发血管壁组织产生光声信号，安置在探头顶端的超声探测器采集来自各个方向的超声回波并传送至计算机，最后重建出血管横截面的二维光吸收系数分布或者初始光声压强分布图像。

图像重建是IVPA成像必不可少的组成部分，其实质就是根据超声探测器采集到的生物组织光声信号求解组织的光吸收系数分布。不同类型的超声探测器和扫描方式对应不同的图像重建算法，IVPA成像的扫描孔径在血管腔内是封闭的，使采集光声信号的方式受到高度限制。目前对IVPA成像导管的研究主要限于单阵元探测器的二维圆周扫描和环形阵列探测器的单方向采集方式。实际上，这两种采集信号的方式是等效的，差别仅在于环形阵列探测器不需要像单阵元探测器一样进行圆周旋转，就可同时接收全方位的光声信号，因此这两种信号采集方式的图像重建算法可以通用。目前，IVPA图像的重建多采用滤波反投影(filtered back-projection，FBP)算法，该算法的优点是原理简单，运算速度快，但存在成像精度低的缺点，因此重建效果不够理想，有必要研究新的重建方法。

技术实现要素：

本发明的目的在于针对现有技术之弊端，提供一种血管内光声图像的时间反演重建方法，以提高成像精度，改善血管内光声图像的重建效果。

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

一种血管内光声图像的时间反演重建方法，所述方法首先以成像导管中心为图像中心建立初始图像；然后将扫描轨迹圆上的各个探测器测量位置看作是时变的点光声信号源，建立超声波在均匀无损的生物组织中的反向传播模型；最后根据建立的超声波反向传播模型，重建出血管横截面的初始光声压强分布图像，所述方法包括以下步骤：

a.建立初始图像：

IVPA的成像平面通过超声探测器并垂直于成像导管，超声探测器的扫描轨迹为位于成像平面内且半径等于导管半径的圆形轨迹，图像重建区域位于扫描轨迹圆的外部；

初始图像A的宽度和高度均为l(单位：mm)，A由M×M个正方形网格组成，相邻网格的间距为Δx＝l/M，各网格点的初始光声压强值均为0，成像平面所在的坐标系为二维笛卡尔直角坐标系XOY，其中坐标原点O是成像导管中心；

b.建立超声波的反向传播模型

超声探测器在测量位置rs处、时刻t∈[0,T]记录的光声压测量值为p'(rs,t)，将扫描轨迹圆上的各个探测器测量位置看作是时变的点光声信号源，用表示成像区域Ω的计算机模拟，在内模拟超声波的传播过程，对Ω内的初始光声压分布进行近似重建，其中，位于成像导管外部，的边界根据时间界限T确定，则超声波在均匀无损的生物组织中的反向传播模型如下：

初始条件为：

式中，c为声速；Δt是时间步长，CFL是仿真精确度和计算速度之间的折中系数；p(r,t)为声压场中位置r∈Ω处在时刻t的声压；p(re,t)是内的非探测器测量位置re处在时刻t的光声压强；p'(rs,T)是超声探测器在测量位置rs处、时刻T记录的光声压测量值；p(rs,t)是声压场中的位置rs处、时刻t的声压；uξ(r,t)为介质在声压场中的位置r处、时刻t的振动速度在成像平面直角坐标系的X轴和Y轴方向上的分量；u(r,t)为介质在声压场中的位置r处、时刻t的振动速度；ρ(r,t)是在声压场中的位置r处、时刻t的声学密度；ρ0是介质密度；ρξ(r,t)为声压场中的位置r处、时刻t的声学密度在成像平面直角坐标系的X轴和Y轴方向上的分量；i是虚数单位；F{...}与F^-1{...}分别是二维傅立叶变换和逆傅立叶变换；κ＝sinc(ckΔt/2)是k空间算子；kξ为在ξ＝(x,y)方向上的空间波数分量；

c.重建血管横截面的初始光声压强分布图像

以Δt为时间步长从t＝0时刻开始依次进行迭代，计算并记录每个测量位置发出的光声波在各个网格点产生的光声压强值，对于某一个网格点来说，它在某一时刻的光声压强值等于该时刻所有测量位置在该网格点产生的光声压强值之和，以t＝T时刻为截止条件，计算此时刻成像平面上每个网格点的光声压强值，即可得到对成像区域Ω内的初始光声压分布的近似重建，最后，将各个网格点的光声压强值转换为灰度矩阵，即可得到血管横截面的灰阶图像。

上述血管内光声图像的时间反演重建方法，在计算每个测量位置发出的光声波在各个网格点产生的光声压强值时，导管内部的网格点不做计算，其光声压强值始终为0。

本发明利用超声探测器采集到的来自于血管壁组织的光声信号数据，在时域内模拟光声信号的反向传播过程，反演得到血管横截面的二维灰阶光声压强分布图像，显示血管壁内部的组织结构。该方法不受公理化推导公式的限制，约束条件少，鲁棒性强，所依赖的假定或初始条件少，而且不易受图像伪影的影响，因此具有更高的成像精度，可获得比较理想的重建效果。仿真实验结果表明，在相同测量位置数下，采用本发明方法重建出的图像的结构相似性指标(structural similarity，SSIM)值可比滤波反投影算法提高约65％。

附图说明

图1是IVPA成像和图像重建示意图，其中图1(a)是IVPA成像示意图；图1(b)是IVPA图像重建示意图；

图2是仿真血管横截面模型；

图3是测量位置数为360时的图像重建结果示意图，其中图3(a)是采用FBP算法的重建图像；图3(b)是采用本发明方法的重建图像。

文中各符号为：A、初始图像；l、图像A的宽度/高度(单位：mm)；M、图像A的宽度/高度对应的离散网格个数；Δx、相邻网格的间距；XOY、成像平面所在的二维直角坐标系；O、坐标系XOY的坐标原点；p(r,t)、声压场中的位置r∈Ω处在时刻t的声压；Ω、成像区域，即光声信号传播的物理区域；uξ(r,t)为介质在声压场中的位置r处、时刻t的振动速度在成像平面直角坐标系的X轴和Y轴方向上的分量；u(r,t)、介质在声压场中的位置r处、时刻t的振动速度；c、声速；ρ(r,t)、声压场中的位置r处、时刻t的声学密度；ρ0、介质密度；p'(rs,t)、超声探测器在测量位置rs处、时刻t∈[0,T]记录的光声压测量值；T、时间界限；成像区域Ω的计算机模拟；re、内的非探测器测量位置；p(re,t)、内的非探测器测量位置re处在时刻t的光声压强；p'(rs,T)、超声探测器在测量位置rs处、时刻T记录的光声压测量值；p(rs,t)、声压场中的位置rs处、时刻t的声压；i、虚数单位；kξ、在ξ＝(x,y)方向上的空间波数分量；F{...}、F^-1{...}、二维傅立叶变换和逆傅立叶变换；Δt、时间步长；κ、k空间算子；ρξ(r,t)、声压场中的位置r处、时刻t的声学密度在成像平面直角坐标系的X轴和Y轴方向上的分量；CFL、仿真精确度和计算速度之间的折中系数。

具体实施方式

下面结合附图对本发明作进一步详述。

本发明首先以成像导管中心为图像中心建立初始图像；然后，建立超声波在均匀无损的生物组织中的反向传播模型；最后，重建出血管横截面的初始光声压强分布图像。具体步骤如下：

(1)建立初始图像

如附图1(a)所示，IVPA的成像平面通过超声探测器并垂直于成像导管。为了简化问题，本发明方法将超声探测器看作理想的点探测器，其扫描轨迹为位于成像平面内且半径等于导管半径的圆形轨迹，图像重建区域位于扫描轨迹圆的外部。

如附图1(b)所示，初始图像A的宽度和高度均为l(单位：mm)，A由M×M个正方形网格组成，相邻网格的间距为Δx＝l/M，各网格点的初始光声压强值均为0。成像平面所在的坐标系为二维笛卡尔直角坐标系XOY，其中坐标原点O是成像导管中心。

(2)建立超声波的反向传播模型

IVPA图像的重建等同于在成像区域内求解t＝0时刻的初始光声压分布。光声信号的本质是超声波，在超声均匀无损传播介质中，超声波传播的物理模型可由以下三个耦合超声方程表示：

式中，p(r,t)为声压场中位置r∈Ω处在时刻t的声压；Ω为成像区域，即光声信号传播的物理区域；u(r,t)为介质在声压场中的位置r处、时刻t的振动速度；c为声速；ρ(r,t)是在声压场中的位置r处、时刻t的声学密度；ρ0是介质密度。

如附图1(b)所示，超声探测器在测量位置rs处、时刻t∈[0,T]记录的光声压测量值为p'(rs,t)，将扫描轨迹圆上的各个探测器测量位置看作是时变的点光声信号源。用表示Ω的计算机模拟，在内模拟超声波的传播过程，对Ω内的初始光声压分布进行近似重建。其中，位于成像导管外部，根据时间界限T可以确定的边界。T值是根据重建区域的大小确定的，即

在图像重建过程中，式(1)的初始条件为

式中，p(re,t)是内的非探测器测量位置re处在时刻t的光声压强；p'(rs,T)是超声探测器在测量位置rs处、时刻T记录的光声压测量值；p(rs,t)是声压场中的位置rs处、时刻t的声压。

本发明采用PSM(pseudospectral method)(Treeby B E,Zhang E,Cox B T.Photoacoustic tomography in absorbing acoustic media using time reversal.Inverse Problems,2010,11(26):115003-115020.)和k空间法(Tabei M,Mast T D,Waag R C.Ak-space method for couple first-order acoustic propagation equations.Journal of Acoustic Society of America,2002,111(1):53-63.)对式(1)中的三个耦合超声方程进行离散化处理，得到：

式中，i是虚数单位；kξ为在ξ＝(x,y)方向上的空间波数分量；F{...}与F^-1{...}分别是二维傅立叶变换和逆傅立叶变换；uξ(r,t)为介质在声压场中的位置r处、时刻t的振动速度在成像平面直角坐标系的X轴和Y轴方向上的分量；Δt是时间步长；κ＝sinc(ckΔt/2)是k空间算子；ρξ(r,t)为声压场中的位置r处、时刻t的声学密度在成像平面直角坐标系的X轴和Y轴方向上的分量；Δt是时间步长，采用k空间法计算得到

式中，Δx是相邻网格的间距，CFL是仿真精确度和计算速度之间的折中系数。

结合式(1)可知，初始声学密度分布为

在图像重建过程中，当r＝re时，即对于内的非探测器测量位置，式(5)中的p(r,t)|t＝0＝p(re,t)|t＝0＝0；当r＝rs时，即对于内的探测器测量位置，式(5)中的p(r,t)|t＝0＝p(rs,t)|t＝0＝p′(rs,T)。

(3)重建血管横截面的初始光声压强分布图像

将式(2)和式(5)作为式(3)的初始条件，以式(4)求得的Δt值为时间步长从t＝0时刻开始依次进行迭代，计算并记录每个测量位置发出的光声波在各个网格点产生的光声压强值(导管内部的网格点不做计算，其光声压强值始终为0)。对于某一个网格点来说，它在某一时刻的光声压强值等于该时刻所有测量位置在该网格点产生的光声压强值之和。以t＝T时刻为截止条件，计算此时刻成像平面上每个网格点的光声压强值，即可得到对成像区域Ω内的初始光声压分布的近似重建。最后，将各个网格点的光声压强值转换为灰度矩阵，即可得到血管横截面的灰阶图像。

分别采用滤波反投影算法和本发明方法对附图2所示的仿真血管横截面模型进行图像重建，结果如附图3(a)和(b)所示，其中成像位置数均为360，图像的结构相似度(Structural Similarity,SSIM)指标分别是0.5049和0.8642，表明本发明方法重建出的图像更接近附图2所示的原始图像。