一种生物磁光声联合内窥成像方法与流程

本发明涉及一种对生物腔体组织进行磁光声联合内窥成像的方法，属于医学成像技术领域。

背景技术：

生物光声(photoacoustic,pa)成像是以生物组织的光声效应为物理基础，即组织吸收短脉冲激光进而发热膨胀产生超声波(即光声信号)。感应式磁声(magnetoacoustictomographywithmagneticinduction,mat-mi)成像以生物组织的磁声效应为物理基础，即将目标体置于静磁场中，并在与静磁场的相同方向上施加磁脉冲激励，目标体中感应出的涡电流与静磁场作用产生洛伦兹力，带电粒子在洛伦兹力的作用下产生振动，从而发射出超声波，再经过超声探头接收后成像。单一的成像技术不能全面、详尽地描述生物组织的结构和功能信息，而超声(ultrasonic,us)、pa和mat-mi成像都是以超声波为载体的声学成像技术，且具有互补的特点，可将它们结合起来，进行磁光声(magneto-photo-acoustic,mpa)联合成像。mpa成像综合了光声信号发射阶段超声检测较高的分辨率以及磁声信号检测时较高的分辨率和灵敏度，可对早期病变组织进行精准的定位和功能成分成像。

目前mpa联合成像都是体外成像，即将超声换能器置于生物体四周，在体外接收超声信号。相对于内窥成像，这种方式不能及时有效地对生物腔体组织(如消化道、肠道和血管等)进行观察和诊断。因此研究一种磁光声联合内窥(endoscopicmagneto-photo-acoustic,empa)成像方法是十分必要的。

技术实现要素：

本发明的目的在于针对现有技术之弊端，提供一种生物磁光声联合内窥成像方法，以获取高分辨率和灵敏度的生物腔体组织图像，帮助医务人员及时有效地对生物腔体组织进行观察和诊断。

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

一种生物磁光声联合内窥成像方法，所述方法在建立多层腔体组织横截面模型的基础上，对腔道内的超声回波成像、光声成像和感应式磁声成像的过程进行数值仿真，得到腔道横截面上组织反射的超声回波信号以及组织产生的光声和磁声信号，然后对三种超声信号进行最优加权求和，得到融合后的联合成像信号。

上述生物磁光声联合内窥成像方法，所述方法包括以下步骤：

a.建立多层腔体组织横截面模型

成像导管位于多层腔体组织横截面模型的中心，超声换能器位于成像导管顶端，模型所在的坐标系是θ-l极坐标系，其中坐标原点是成像导管顶端中心，θ是极角，l是极径，水平向右的方向为l轴正方向，忽略超声换能器的孔径效应，将超声换能器看作理想的点探测器，其扫描轨迹为平行于成像平面的圆形轨迹，以模型的中心为起始点，将模型等角度划分为n份，每一份近似为层与层之间平行的多层腔体组织，导管处的成像角度为

θi＝360(i-1)/n

其中i＝1,2,…,n，θi对应的成像区域的角度范围为[θia,θib]，其中θia＝θi-180/n，θib＝θi+180/n；

b.对成像组织反射的超声回波信号以及组织产生的光声信号和磁声信号进行数值仿真；

①超声回波信号：

从腔体中心沿径向发射超声脉冲，根据不同成分组织的声阻抗差异值以及超声探测器的冲激响应，仿真得到组织反射/散射的超声回波信号；

②光声信号：

从腔体中心对周围组织沿径向发射激光脉冲，组织由于光声效应产生光声信号，根据不同成分组织的光吸收系数和散射系数，结合蒙特卡罗模拟和光声波动方程，仿真得到组织产生的光声信号；

③磁声信号：

沿腔体的轴向施加静磁场和脉冲磁激励，组织由于磁声效应产生磁声信号，根据不同成分组织的电导率，结合磁声波动方程，仿真得到组织产生的磁声信号；

c.对超声回波信号、光声信号和磁声信号进行融合

对于成像组织中的角度θi(i＝1,2,…,n)、位置r处的超声回波信号光声信号和磁声信号采用下式进行融合：

其中为超声回波、光声和磁声信号的融合信号，即联合成像信号，w1i、w2i、w3i分别为超声回波、光声和磁声信号的加权因子，且有w1i+w2i+w3i＝1。

上述生物磁光声联合内窥成像方法，为了使超声回波、光声和磁声信号的融合信号的均方误差最小，超声回波、光声和磁声信号的加权因子w1i、w2i、w3i应取最优值最优值的计算公式为：

其中k为超声回波、光声和磁声信号的测量次数，为总均方误差最小时超声回波、光声和磁声信号的最优加权因子，为的方差，为的方差，为的方差，ruu为的自相关系数，rup为和的互相关系数，rpp为的自相关系数，rpm为和的互相关系数，rmm为的自相关系数，rmu为和的互相关系数。

本发明在信号层对超声换能器分时接收的腔体组织反射/散射的超声回波信号以及组织产生的光声信号和磁声信号进行融合，同传统方法相比，本方法得到的联合成像信号可较多地保留组织的形态结构和成分信息，据此重建出的组合图像具有极高的空间分辨率、对比度、灵敏度和对比分辨率，能准确显示腔道壁内各组织的位置、形态及其功能成分。

附图说明

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

图1是含有钙化斑块的血管横截面模型示例；

图2是empa成像导管在角度θi处分时接收超声回波信号、光声信号和磁声信号的示意图；

图3是将图2中角度θi处对应的成像区域近似为多层血管壁组织的示意图。

文中各符号为：θ、l为θ-l平面极坐标系的极角和极径，其中成像导管位于坐标原点(腔体中心)，水平向右的方向为l轴正方向；n为腔体横截面模型被等角度分割的份数；θi为第i个成像角度；θia、θib为θi对应的成像区域角度范围的下限和上限，其中i＝1,2,…,n；f(r,θ)为超声探测器在位置r、角度θ处采集的超声回波信号；t(r,θ)为成像组织在位置r、角度θ处的声阻抗差异函数；h(r,θ)为超声探测器的点扩散函数，即冲激响应；σr为发射超声波的脉冲宽度；σθ为发射超声波束的宽度；k0为波数；f0为发射超声波的中心频率；g(r,θ)为均值为0、标准差为1的高斯白噪声；e(r,θ)为位置r、角度θ处组织的声阻抗差异值；为角度θi、位置r处采集的超声回波信号，其中i＝1,2,…,n；为哈密顿算子；p(r,t)为时刻t、位置r处的声压；c为生物组织中的声速；a(r)为光能量沉积分布函数；β为组织的体积膨胀温度系数；cp为组织的比热容；i(t)为光声成像中入射激光脉冲的强度；(j,k)为组织上的点r在极坐标系θ-l中的坐标；δθ、δl为θ轴和l轴上的单位长度；δt为离散时间间距；n为离散时刻；pⁿ(j,k)为时刻n、位置(j,k)处的声压；vθⁿ(j,k)和vlⁿ(j,k)为时刻n、位置(j,k)处的质点分别在θ方向和l方向的振动速度；ρ为组织密度；c(j,k)为位置(j,k)处的声速；a(j,k)为位置(j,k)处的光能量沉积值；iⁿ为时刻n的激光脉冲强度；cmax为组织中声速的最大值；λmin为超声波的最小波长；为角度θi、位置r处的光声信号，其中i＝1,2,…,n；j为感应电流密度；b0为静磁场强度；j(j,k)为位置(j,k)处的感应电流密度；b0(j,k)为位置(j,k)处的静磁场强度；为角度θi、位置r处的磁声信号，其中i＝1,2,…,n；k为超声回波、光声和磁声信号的测量次数；ruu为的自相关系数；rpp为的自相关系数；rmm为的自相关系数；rup为和的互相关系数；rpm为和的互相关系数；rmu为和的互相关系数；为的方差；为的方差；为的方差；w1i、w2i、w3i为超声回波、光声和磁声信号的加权因子；为超声回波、光声和磁声信号的融合信号；为融合信号的总均方误差；为总均方误差的最小值；为总均方误差最小时超声回波、光声和磁声信号的最优加权因子。

具体实施方式

磁光声联合内窥(endoscopicmagneto-photo-acoustic,empa)成像是用同一个成像系统在生物腔道内部同时进行超声、光声和感应式磁声成像。超声换能器直接在腔道内采集组织反射、散射或者产生的超声波信号，再经计算机对其进行融合后获得组合图像。

一、建立多层腔体组织横截面模型：

如附图1所示，以血管横截面模型为例，成像导管位于模型的中心，超声换能器位于成像导管顶端，环绕导管由内向外沿径向依次是血管内腔、粥样硬化斑块(钙化、脂质、纤维或混合斑块)、血管壁内膜/中膜和外膜。模型所在的坐标系是θ-l极坐标系，其中坐标原点是成像导管中心，θ是极角，l是极径，水平向右的方向为l轴正方向。本发明方法忽略超声换能器的孔径效应，将超声探测器(即超声换能器)看作理想的点探测器，其扫描轨迹为平行于成像平面的圆形轨迹。

如附图2所示，以血管模型的中心为起始点，将模型等角度划分为n份，每一份近似为层与层之间平行的多层血管壁组织(如附图3所示)。对模型分别施加超声脉冲、激光脉冲和激励磁场(如附图2所示，包含静磁场与脉冲磁场，其中静磁场的方向与腔体的轴向重合)。导管处的成像角度为

θi＝360(i-1)/n(1)

其中i＝1,2,…,n。θi对应的成像区域的角度范围为[θia,θib]，其中θia＝θi-180/n，θib＝θi+180/n。

二、仿真腔体组织的超声回波信号：

超声内窥成像是采用超声束在腔道内周向旋转扫描，并通过超声脉冲反射法检测病变组织。生物软组织是有层次的，不同层次的组织成分不同，表现为声学特征参量的差异，超声波在组织中传播时，当遇到不同组织之间的分界面时，会产生反射回波，其中包含了不同组织的位置和结构信息。

超声探头从腔体横截面模型的中心沿径向发射超声信号，为简化问题，不考虑组织的吸收衰减和高阶回波，则：

f(r,θ)＝h(r,θ)*t(r,θ)(2)

其中，f(r,θ)是超声探测器在位置r、角度θ处采集的超声回波信号；h(r,θ)是是超声探测器的点扩散函数，即冲激响应，它是时不变且可分离的，表示为

式中，σr是发射超声波的脉冲宽度，σθ是发射超声波束的宽度，波数k0＝2πf0/c，c是组织中的声速，f0是发射超声波的中心频率。式(2)中的t(r,θ)是成像组织在位置r、角度θ处的声阻抗差异函数：

t(r,θ)＝g(r,θ)*e(r,θ)(4)

其中，g(r,θ)是均值为0、标准差为1的高斯白噪声；e(r,θ)是位置r、角度θ处组织的声阻抗差异值，若组织的e值较大，则表明其声阻抗差异较大，产生的回波信号也较强。

由式(2)得到超声探测器在n个角度采集的超声回波信号，将角度θi、位置r处采集的超声回波信号记为(i＝1,2,…,n)。

三、仿真短脉冲激光作用于腔体组织产生的光声信号：

从腔体横截面的中心沿径向发射短脉冲激光，根据组织的光吸收系数和散射系数，通过组织内光的蒙特卡罗模拟，获得光能量沉积分布函数。由光声效应导致组织产生的光声信号(其实质是超声波)满足光声波动方程

其中，是哈密顿算子；p(r,t)是时刻t、位置r处的声压；c是生物组织中的声速；a(r)是位置r处的光能量沉积分布；β是组织的体积膨胀温度系数，cp是组织的比热容；i(t)是入射激光脉冲的强度。

结合超声波的声压、质点振动速度和组织密度三个物理量之间的关系，采用时域有限差分法对式(5)进行离散化处理得到

其中，(j,k)是组织上的点r在极坐标系θ-l中的坐标；δθ和δl分别是θ轴和l轴的单位长度；δt是离散时间间距；n是离散时刻；pⁿ(j,k)是时刻n、位置(j,k)处的声压；和是时刻n、位置(j,k)处的质点分别在θ方向和l方向的振动速度；ρ是组织密度；c(j,k)是位置(j,k)处的声速；a(j,k)是位置(j,k)处的光能量沉积值；iⁿ是时刻n的激光脉冲强度。

式(6)需满足courant稳定性条件：

其中，cmax是组织中声速的最大值；λmin是超声波的最小波长。

由式(6)可以得到n个角度的光声信号，将角度θi、位置r处的光声信号记为(i＝1,2,…,n)。

四、仿真静磁场和脉冲磁场作用于腔体组织产生的磁声信号：

在mat-mi成像中，生物组织因受到外界脉冲磁场的激励而产生感应涡流，感应涡流与外界施加的静磁场作用产生洛伦兹力，组织中的带电粒子在洛伦兹力的作用下产生机械振动，并以超声波的形式向外传播，即磁声信号。

为简化起见，将组织近似为理想的均匀流体，忽略其黏性和不可压缩性，其初始状态是静止的，且忽略热交换，假设声波传播过程为绝热过程。基于以上假定，组织产生的磁声信号满足磁声波动方程

其中，p(r,t)是时刻t、位置r处的声压；c是生物组织中的声速；j是感应电流密度，b0是静磁场强度。

式(8)的时域有限差分形式为：

其中，(j,k)是组织上的点r在极坐标系θ-l中的坐标；δθ和δl分别是θ轴和l轴的单位长度；δt是离散时间间距；n是离散时刻；pⁿ(j,k)是时刻n、位置(j,k)处的声压；和是时刻n、位置(j,k)处的质点分别在θ方向和l方向的振动速度；ρ是组织密度；c(j,k)是位置(j,k)处的声速；j(j,k)是位置(j,k)处的感应电流密度；b0(j,k)是位置(j,k)处的静磁场强度。

式(9)也需满足式(7)中的courant稳定性条件。由式(9)可以得到n个角度的磁声信号，将角度θi、位置r处的磁声信号记为(i＝1,2,…,n)。

五、融合超声回波、光声信号和磁声信号：

本发明方法在总均方误差最小的条件下，对于成像组织在角度θi(i＝1,2,…,n)、位置r处的超声回波信号、光声信号和磁声信号的测量值，以自适应的方式寻找各信号对应的最优加权因子，得到最优的融合结果。具体步骤如下：

设和的测量次数均为k，测量值的方差分别是和具体计算方法如下：

对于和的自相关系数为

和的互相关系数为

则

对于和的自相关系数为

其中

和的互相关系数为

则

对于和的自相关系数为

和的互相关系数为

则

超声回波、光声和磁声信号的加权因子分别为w1i、w2i和w3i，且满足

w1i+w2i+w3i＝1(19)

融合后的信号为

总均方误差为

通过对式(21)中的多元二次函数求极值，得到总均方误差的最小值

总均方误差最小时，超声回波、光声和磁声信号的最优加权因子分别是：

将式(23)代入式(20)即可得到融合后的信号。该融合算法无需有关待融合信号的任何先验知识，只需根据待融合信号的测量值，估计各信号方差的变化，及时调整参与融合的各信号的加权因子，得到均方误差最小的融合信号。融合后信号的均方误差不仅小于单个信号的均方误差，而且融合结果在精度、容错性方面均优于传统的平均值估计算法。