一种磁感应磁声内窥图像的建模与仿真方法

摘要

一种磁感应磁声内窥图像的建模与仿真方法，所述方法首先建立生物腔体组织的横截面模型和生物腔体组织横截面的电磁特性参数模型；然后应用有限元分析的方法，对腔体组织在磁激励场中产生感应电流的过程进行仿真，再由感应电流仿真腔体组织在静磁场中产生声源的过程；之后根据声源的分布仿真出腔体组织产生的磁声信号；最后根据磁声信号重建出生物腔体的横截面图像。本发明可灵活地调整生物腔体组织横截面模型，改变模型中腔体组织的变异类型及程度，设置多层腔体壁组织的电磁特性参数和声学参数，精确地仿真声速不均匀的多层生物腔体组织的声场，因此可以为成像算法和图像后处理算法的研究和性能测试等提供足够的样本图像数据。

说明书

技术领域

本发明涉及一种医学成像技术，特别是对生物腔体的磁感应磁声内窥横截面图像进行建模与仿真。

背景技术

生物组织磁感应磁声内窥(Endoscopic Magnetoacoustic Tomography with Magnetic Induction,EMAT-MI)成像属于功能成像，其成像的物理基础是生物组织的磁声效应，即生物组织变异部位与非变异部位的电导率不同，导致产生的感应涡流不同，从而在静磁场中产生不同强度的超声信号(即磁声信号)。EMAT-MI结合了生物电阻抗成像的高对比度和超声扫描成像的高空间分辨率，可得到生物腔体组织(含变异组织)内部的形态结构、组织类型和成分等信息。然而EMAT-MI成像技术尚处于实验室研究阶段，未大规模应用于实际，而改善成像仪器的结构和校准参数、优化图像分析和处理算法等都是建立在分析大量样本图像数据的基础之上的。因此，为了获取足够多的样本图像数据，有必要寻找一种对生物腔体组织磁感应磁声内窥图像进行建模与仿真的方法。

发明内容

本发明的目的在于针对现有技术之弊端，提供一种磁感应磁声内窥图像的建模与仿真方法，为成像和图像后处理的研究和性能测试等提供数据源。

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

一种磁感应磁声内窥图像的建模与仿真方法，所述方法首先建立生物腔体组织的横截面模型和生物腔体组织横截面的电磁特性参数模型；然后应用有限元分析的方法，对腔体组织在磁激励场中产生感应电流的过程进行仿真，再由感应电流仿真腔体组织在静磁场中产生声源的过程；之后根据声源的分布仿真出腔体组织产生的磁声信号；最后根据磁声信号重建出生物腔体的横截面图像。

上述磁感应磁声内窥图像的建模与仿真方法，具体处理按以下步骤进行：

a.建立生物腔体组织的横截面模型：

成像导管位于生物腔体组织横截面模型的中心，接收磁声信号的超声探测器位于成像导管顶端，将超声探测器看作理想的点换能器，其扫描轨迹为平行于成像平面、半径趋近于0的圆形轨迹；

b.建立生物腔体组织横截面的电磁特性参数模型：

以生物腔体组织横截面模型的中心为起始点，将模型等角度划分为m份，每一份近似为多层组织，对模型施加脉冲磁激励，并对每一份组织接收到的磁声信号进行仿真，成像导管上的超声探测器所处的角度为

θ_i＝360(i-1)/m

其中，i＝1,2,...,m，对应的成像区域的角度范围为[θ_ia,θ_ib]，其中θ_ia＝θ_i-180/m，θ_ib＝θ_i+180/m。

确定每个成像区域的每层组织的电导率和厚度参数，形成生物腔体组织横截面的电磁特性参数模型。以血管的横截面模型为例，多层血管壁组织的参数设置示例如表1所列。

c.仿真EMAT-MI成像中声源产生的过程：

①应用仿真软件构造亥姆霍兹(Helmholtz)线圈，并对亥姆霍兹线圈施加高斯脉冲电流作为激励源，产生交变磁场B₁，沿亥姆霍兹线圈轴向施加一个稳恒匀强磁场B₀，多层生物腔体组织的横截面模型同轴放置于亥姆霍兹线圈中间；稳恒磁场和亥姆霍兹线圈的参数设置示例如表2所列。

②以空气为背景域，以自由四面体为单元对亥姆霍兹线圈以及生物腔体组织的横截面模型进行划分，进而采用有限元分析的方法仿真得到生物腔体组织中感应涡流J的分布；

③根据感应涡流J仿真出生物腔体组织的声源的分布，其中，“×”表示向量积；

d.根据声源的分布仿真生物腔体组织产生的磁声信号，得到超声探测器在时刻t、角度θ_i、位置r处接收到的多层生物腔体组织产生的磁声信号的声压p_i(r,t)(i＝1,2,...,m)；

e.重建极坐标系下的EMAT-MI图像：

①根据步骤d得到的m个磁声信号的声压p_i(r,t)(i＝1,2,...,m)计算出

$\frac{{\partial }^2{p}_i(r,t)}{\partial t^2}{|}_{t=|r_i-r|/c_s},$

式中，c_s是超声波的波速，r_i是θ-l平面中超声探测器所在位置(与成像导管在X-Y平面中的成像角度θ_i相对应)；

②由下式重建出在角度θ_i处的声源分布：

$▿·{(j\times B_0)}_i\approx -\frac{1}{2{\mathrm{\pi c}}_s^3}\frac{q·(r_i-r)}{|r_i-r{|}^2}\frac{{\partial }^2{p}_i(r,t)}{\partial t^2}{|}_{t=|r_i-r|/c_s};$

式中，q是r_i处的单位矢量；

③将上式代入下式，得到腔体组织在角度θ_i处的电导率分布：

${\sigma }_i\left(r\right)\approx -\frac{▿·{(J\times B_0)}_i}{B_0·B_1},$

将σ_i(r)作为角度θ_i、位置r处的腔体组织横截面极坐标视图的灰度值，从而得到极坐标下的EMAT-MI图像；

f.图像的坐标转换：

将极坐标下的EMAT-MI图像转换为X-Y直角坐标系下的横截面视图，具体方法如下：

建立X-Y平面直角坐标系，坐标原点是成像导管中心，水平向右的方向为X轴正方向，垂直于X轴向上的方向为Y轴正方向，设极坐标系中的一点(j,k)的灰度值为f(j,k)，该点在X-Y坐标系中的对应点的坐标为(j′,k′)，灰度值为g(j′,k′)，其中j∈[0,2π]，k∈[0,d]，j′∈[-d,d]，k′∈[-d,d]，d为极坐标视图中极径的最大值，则有：

g(j′,k′)＝f(j,k)，

其中

$\left(\begin{array}{c}{j}^{\prime }=kcos\left(j\right)\\ k^{\prime }=ksin\left(j\right)\end{array}\right).$

本发明可灵活地调整生物腔体组织横截面模型，改变模型中腔体组织的变异类型及程度，设置多层腔体壁组织的电磁特性参数和声学参数，精确地仿真声速不均匀的多层生物腔体组织的声场，因此可以为成像算法和图像后处理算法的研究和性能测试等提供足够的样本图像数据。

附图说明

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

图2是EMAT-MI成像导管上的超声探测器在角度θ_i处接收磁声信号的示意图；

图3是将图1中的血管横截面模型等角度划分后，将其中的一份近似为多层组织的示意图；

图4是θ-l坐标系中的多层生物腔体组织的示意图；

图5是亥姆霍兹(Helmholtz)线圈激励示意图，图中R是线圈半径。

表1是多层血管壁组织参数设置举例。

表2是稳恒磁场和亥姆霍兹线圈的参数设置示例。

文中各符号为：X、Y、腔体横截面模型所在的X-Y平面直角坐标系的横轴和纵轴，其中，成像导管的中心为坐标原点，水平向右的方向为X轴正方向，垂直于X轴向上的方向为Y轴正方向；m、横截面模型被等角度分割的总份数；θ_i、成像导管的第i个成像角度，其中，i＝1,2,...,m；θ_ia、θ_ib、成像导管在角度θ_i处进行成像时对应的角度范围的上、下限；θ、l、θ-l坐标系的横轴和纵轴；B₁、激励磁场的磁通密度；B₀、稳恒磁场的磁通密度；J、生物腔体组织中的感应涡流；×、向量积；哈密顿算子；j、k、θ-l坐标系沿θ方向和l方向的单位向量；t、时间；r、θ-l坐标系中的一点；p_i(r,t)、超声探测器在时刻t角度θ_i位置r处接收到的多层生物腔体组织产生的磁声信号的声压，其中i＝1,2,...,m；J_i(r)、生物腔体组织在角度θ_i位置r处产生的感应涡流，其中i＝1,2,...,m；B_0i(r)、生物组织在角度θ_i位置r处的激励磁场和稳恒磁场的叠加磁通密度，其中i＝1,2,...,m；η(t)、随时间变化的激励电流信号；c_s、超声波的波速；ρ、生物组织的密度；v_li(r)和v_θi(r)分别是角度θ_i方向位置r处的质点在l方向和θ方向的振动速度，其中i＝1,2,...,m；(j,k)、位于θ-l平面内腔体组织上一点的坐标；Δθ、Δl、坐标轴θ和坐标轴l方向的的离散空间间距；Δt、离散时间间距；n、离散时间点；pⁿ(j,k)、位置(j,k)的质点在时刻n产生的磁声信号的声压；时刻n位置(j,k)的质点在θ方向和l方向的振动速度；ρ(j,k)、生物组织在点(j,k)处的密度；c_s(j,k)、磁声信号在位置(j,k)的传播速度；r_i、θ-l平面中超声探测器所在位置(与成像导管在X-Y平面中的成像角度θ_i相对应)；q、超声探测器所在位置r_i处的单位矢量；σ_i(r)、生物腔体组织在角度θ_i位置r处的电导率，其中i＝1,2,...,m；(j′,k′)、θ-l坐标系中的一点(j,k)在X-Y直角坐标系中对应点的坐标；f(j,k)、g(j′,k′)、点(j,k)和点(j′,k′)的灰度值；d、极坐标视图极径的最大值。

具体实施方式

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

本发明方法的步骤包括：

(1)建立多层生物腔体组织横截面的形态模型：

以血管横截面形态模型为例，如附图1所示，模型包括成像导管(接收磁声信号的超声探测器位于成像导管顶端)、管腔、血管壁内膜/中膜、外膜和斑块(脂质斑块、纤维斑块、钙化斑块或混合斑块)五部分。其中成像导管位于模型中心，周围由内向外依次为管腔、斑块、血管壁内膜/中膜和外膜。可根据组织变异的类型(如钙化、纤维化或者脂质斑块)和大小，以及血管内腔、血管壁内膜/中膜、外膜的厚度建立不同的血管横截面形态模型。本发明方法忽略超声探测器的孔径效应，将其看作理想的点换能器，其扫描轨迹为平行于成像平面、半径趋近于0的圆形轨迹。模型所在的坐标系为X-Y平面直角坐标系，其中坐标原点是成像导管中心，水平向右的方向为X轴正方向，垂直于X轴向上的方向为Y轴正方向。

(2)建立多层生物腔体组织横截面的电磁特性参数模型：

首先，如附图2所示，以生物腔体组织横截面形态模型的中心为起始点，将模型等角度划分为m份，对模型施加脉冲磁激励，并对每一份组织接收到的磁声信号进行仿真。成像导管上的超声探测器所处的成像角度为

θ_i＝360(i-1)/m(1)

其中，i＝1,2,...,m。对应的成像区域的角度范围为[θ_ia,θ_ib]，其中θ_ia＝θ_i-180/m，θ_ib＝θ_i+180/m。例如，当m＝360时，将附图1中含有脂质斑块的血管横截面形态模型以X-Y坐标系的原点为中心等分为360份，超声探测器对每一份组织采集磁声信号的过程中所处的角度为θ₁＝0°,θ₂＝1°,...,θ_i＝360(i-1)/360,...,θ₃₆₀＝359°。

然后，如附图4所示，在θ-l坐标系中，每一份多层组织的表面平行于θ轴且垂直于l轴。其中θ轴正方向为水平向右的方向，l轴正方向为垂直于θ轴向上的方向，表示多层组织的厚度。确定每个成像角度对应的组织中每层的电导率和厚度参数，形成多层生物腔体组织的电磁特性参数模型。以血管的横截面模型为例，多层血管壁组织的参数设置示例如表1所列。不同成像角度对应的多层生物腔体组织仅厚度参数不同，图4中的生物腔体组织分为a、b、c、d、e五层，若形态模型中不包含斑块，则c、d层的厚度参数为0。

(3)仿真EMAT-MI成像中声源产生的过程：

本发明采用有限元分析的方法进行声源的仿真，具体步骤如下：

首先，应用COMSOL软件构造亥姆霍兹(Helmholtz)线圈，如附图5所示，并施加高斯脉冲电流作为激励源产生交变磁场B₁，沿线圈轴向施加一个稳恒匀强磁场B₀，多层生物腔体组织模型置于线圈中间，其轴向与亥姆霍兹线圈的轴向同向；然后，以空气为背景域，以自由四面体为单元对亥姆霍兹线圈以及生物腔体组织模型进行划分，进而仿真得到生物腔体组织中感应涡流J的分布；最后，根据感应涡流J仿真出生物腔体组织的声源的分布，其中，“×”表示向量积。稳恒磁场和亥姆霍兹线圈的参数设置示例如表2所列。

表1

表2

(4)仿真生物腔体组织产生的磁声信号：

生物组织产生的磁声信号的本质是超声波，描述磁声信号在声学均匀介质中传播的物理模型为：

$\left(\begin{array}{c}{▿}^2{p}_i(r,t)-\frac{1}{c_s^2}\frac{{\partial }^2{p}_i(r,t)}{\partial t^2}=-▿·[J_i\left(r\right)\times B_{0i}\left(r\right)]\\ ▿=\frac{\partial }{\partial \theta }j+\frac{\partial }{\partial l}k\end{array}\right)---\left(2\right)$

其中，i＝1,2,...,m；为哈密顿算子；“·”为向量的点积；j和k分别为θ-l坐标系的θ方向和l方向的单位向量；t为时间；r为θ-l坐标系中的一点，如附图4所示；c_s为超声波的波速；p_i(r,t)为超声探测器在时刻t角度θ_i位置r处接收到的多层生物腔体组织产生的磁声信号的声压；J_i(r)为生物组织在角度θ_i位置r处产生的感应涡流；B_0i(r)为生物组织在角度θ_i位置r处稳恒磁场的磁通密度。

生物组织的声学特性与流体相似，在理想流体中，声波的传播满足牛顿第二定律、质量守恒定律和物态方程。设随时间变化的激励电流信号为η(t)，则由声波的连续性方程、运动方程和物态方程(张海澜.理论声学.北京：高等教育出版社.2007.)可将式(2)改写成：

$\left(\begin{array}{c}\frac{\partial p_i(r,t)}{\partial t}=-{{\mathrm{\rho c}}_s}^2(\frac{\partial v_{li}\left(r\right)}{\partial l}+\frac{\partial v_{\theta i}\left(r\right)}{\partial \theta })+{c_s}^2\eta \left(t\right)[▿\left(J_i\right(r)\times B_{0i}(r\left)\right)]\\ \frac{\partial v_{\theta i}\left(r\right)}{\partial t}=-\frac{1}{\rho }\frac{\partial p_i(r,t)}{\partial \theta }\\ \frac{\partial v_{li}\left(r\right)}{\partial t}=-\frac{1}{\rho }\frac{\partial p_i(r,t)}{\partial l}\end{array}\right)---\left(3\right)$

其中，i＝1,2,...,m；p_i(r,t)是声波的声压；ρ是生物组织的密度；v_li(r)和v_θi(r)分别是角度θ_i方向位置r处的质点在l方向和θ方向的振动速度。

采用时域有限差分法(Yee K S.Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media.IEEE Trans.Antennas Propag,1966,14(3):302-307.)和Cartesian网格划分法(Sakaguchi T,Hirano T,Watanabe Y,et al.Inner head acoustic field for bone-conducted sound calculated by finite-difference time-domain method.Japanese journal of applied physics,2002,41(5S):3604-3608.)对方程(3)中的各个物理量离散化得到：

$\left(\begin{array}{c}{p}^{n+1}(j,k)=p^n(j,k)-\rho (j,k){c_s}^2(j,k)\Delta t\\ \left\{\frac{1}{\Delta \theta }\right[v_{\theta }^{n+\frac{1}{2}}(j+\frac{1}{2},k)-v_{\theta }^{n+\frac{1}{2}}(j-\frac{1}{2},k)]+\\ \frac{1}{\Delta l}[v_l^{n+\frac{1}{2}}(j,k+\frac{1}{2})-v_l^{n+\frac{1}{2}}(j,k-\frac{1}{2})]\}+\\ {c_s}^2(j,k){\eta }^n[▿\left(J\right(j,k)\times B_0(j,k\left)\right)]\Delta t\\ v_{\theta }^{n+\frac{1}{2}}(j-\frac{1}{2},k)=v_{\theta }^{n-\frac{1}{2}}(j-\frac{1}{2},k)-\frac{\Delta t}{\rho ·\Delta \theta }[p^n(j,k)-p^n(j-1,k)]\\ v_l^{n+\frac{1}{2}}(j,k-\frac{1}{2})=v_l^{n-\frac{1}{2}}(j,k-\frac{1}{2})-\frac{\Delta t}{\rho ·\Delta l}[p^n(j,k)-p^n(j,k-1)]\end{array}\right)---\left(4\right)$

式中，(j,k)表示位于θ-l平面内腔体组织上一点的坐标；Δθ和Δl分别表示坐标轴θ和坐标轴l方向的离散空间间距；Δt表示离散的时间间距；n表示离散时间点；pⁿ(j,k)是位置(j,k)的质点在时刻n产生的磁声信号的声压；和是时刻n位置(j,k)的质点在θ方向和l方向的振动速度；ρ(j,k)是生物组织在点(j,k)处的密度；c_s(j,k)是磁声信号在位置(j,k)的传播速度。

根据式(4)即可仿真得到超声探测器在时刻t、角度θ_i、位置r处接收到的多层生物腔体组织产生的磁声信号的声压p_i(r,t)(i＝1,2,...,m)。

(5)重建极坐标下的EMAT-MI图像：

图像重建的实质是由超声探测器接收到的磁声信号计算出生物组织的电导率分布。具体方法如下：

首先，根据步骤(4)得到的m个磁声信号p_i(r,t)(i＝1,2,...,m)计算出

$\frac{{\partial }^2{p}_i(r,t)}{\partial t^2}{|}_{t=|r_i-r|/c_s}---\left(5\right)$

然后，由下式重建出在角度θ_i处的声源分布：

$▿·{(j\times B_0)}_i\approx -\frac{1}{2{\mathrm{\pi c}}_s^3}\frac{q·(r_i-r)}{|r_i-r{|}^2}\frac{{\partial }^2{p}_i(r,t)}{\partial t^2}{|}_{t=|r_i-r|/c_s}---\left(6\right)$

式中，i＝1,2,...,m；r_i是θ-l平面中超声探测器所在位置(与成像导管在X-Y平面中的成像角度θ_i相对应)；q是r_i处的单位矢量。

最后，将式(6)代入下式重建出腔体组织在角度θ_i处的电导率分布：

${\sigma }_i\left(r\right)\approx -\frac{▿·{(J\times B_0)}_i}{B_0·B_1}---\left(7\right)$

其中，i＝1,2,...,m。将σ_i(r)作为角度θ_i、位置r处的腔体组织横截面极坐标视图的灰度值。

(6)图像的坐标转换：

将步骤(5)中获得的θ-l坐标系中的极坐标视图转换为X-Y直角坐标系下的横截面视图。具体方法如下：

设θ-l坐标系中的一点(j,k)的灰度值为f(j,k)，该点在X-Y坐标系中的对应点的坐标为(j′,k′)，灰度值为g(j′,k′)，其中j∈[0,2π]，k∈[0,d]，j′∈[-d,d]，k′∈[-d,d]，d为极坐标视图纵坐标的最大值。那么有：

g(j′,k′)＝f(j,k)(8)

其中

$\left(\begin{array}{c}{j}^{\prime }=kcos\left(j\right)\\ k^{\prime }=ksin\left(j\right)\end{array}\right)---\left(9\right)$