基于快速体积分方程和磁共振的人体电磁特性反演方法

摘要

基于快速体积分方程和磁共振的人体电磁特性反演方法，涉及核磁共振成像。结合玻恩迭代算法或变分玻恩迭代算法或变形玻恩迭代算法、快速傅里叶变换、稳定双共轭梯度算法以及共轭梯度算法，进行反复的正演和反演迭代至结果收敛，求解人体电磁特性参数，实现人体电磁特性参数磁共振成像的方法，得到人体组织电磁特性参数分布的二维或三维图像，呈现组织内各区域电磁特性参数差异的同时，也可对其进行定量研究，可用于医学上的疾病研究和指导临床诊断、治疗。

说明书

技术领域

本发明涉及核磁共振成像，尤其是涉及一种基于快速体积分方程和磁共振的人体电磁特性反演方法。

背景技术

1991年Haacke第一次提出了基于核磁共振的电特性成像(EPT)。电特性成像与核磁共振一样都是非侵入式的成像方法，它利用核磁共振的数据反演出组织的电特性参数分布。磁共振检测伴随着热效应，这可能给人体带来伤害，尤其是未来要发展的高场磁共振检测，热效应与人体的电特性参数分布有关，了解电特性参数的分布可以提前评估出磁共振检测给人体带来的热效应及其影响，这对于核磁共振成像的发展有很大的意义。癌症由于其高死亡率和低治愈率而受到人们的广泛关注，有研究表明，癌症细胞比正常细胞的电特性参数高很多，比如乳癌细胞比正常细胞大200％，而膀胱癌细胞比正常细胞大100％，利用电特性成像就可以更清晰明了地分辨出肿瘤和正常组织，这有利于癌症的早期发现和及时理疗。目前，基于核磁共振的电特性成像方法已经有多种，比如正交鸟笼线圈法、多通道传输接收法、局部麦克斯韦成像法等等。但是这些方法都有一定的局限性，它们都基于麦克斯韦等式，大多都假设电特性参数的分布是局部均匀的，没有将电磁场边界条件考虑在内，虽然简化了计算，但使得电特性参数在组织的边界处会变得不可靠，且对于噪声也更加敏感。Balidemaj等人虽提出了利用核磁共振系统测量的数据进行迭代反演的对比源反演方法，有效地抑制边界出现的错误，但作者只重建了二维模型，且没有将射频屏蔽考虑在内。

21世纪初定量磁化率成像(QSM)被提出，利用MRI方法获取磁化率定量图像，对磁化率做定量分析，可以为生物医学研究提供一种非创伤性手段，而且对组织内部顺磁性铁含量的测量，有利于对脑部血管疾病和神经病变疾病的诊断和治疗。QSM需要对相位信息进行解缠绕及去除背景场的预处理来获得反映局部磁场变化的场图，再结合重建算法重建出磁化率图像。对于QSM方法，在病态逆问题的求解情况下，背景场去除效果不理想的话会对磁化率的求解产生较大影响，甚至会导致最终的磁化率成像结果不可用，目前两种去除背景场较好的方法是复杂谐波伪影去除法(SHARP)和偶极场投影法(PDF)。由场图信息重建磁化率图像 是一个不适定逆问题，稳定求解和精确定量是关键问题。常用的QSM重建方法有多方向采样磁化率计算方法(COSMOS)、贝叶斯正则化方法、k空间加权微分法(WKD)等。由于成像仪器腔体限制和病人舒适角度问题，很难获得多个不同摆放方向的脑部成像，这极大限制了COSMOS方法的临床应用。贝叶斯正则化方法，随着匹配噪声参数增大，磁化率重建图像清晰度变好，显示出很好的去噪能力。WKD方法成像能显示出脑组织局部细节，且没有出现明显的条状伪影，但成像参数(如场强、回波时间、翻转角等)会对相位处理和磁化率分布图的精度产生影响。

发明内容

本发明的目的在于针对现有技术中使用的电特性成像(EPT)算法在组织边界处值不可靠，对于噪声更加敏感，导致成像分辨率较差，以及磁共振(MR)成像中幅值信息和相位信息分开成像，且利用相位信息的磁化率成像需要解缠绕和去除背景场处理等问题，提供一种基于快速体积分方程和磁共振的人体电磁特性反演方法。

本发明包括以下步骤：

1)测量实验数据，利用磁共振的影像技术测量人体，得到人体的B₁⁺场，所述B₁⁺场是磁共振射频场的正旋场；

2)进行迭代计算，具体步骤为：

①假定初始的电特性对比度χ^E；(0)(r)和磁化率χ^H；(0)(r),利用方程(1)和(2)，进行正演计算，解出E⁽⁰⁾(r)和H⁽⁰⁾(r)；

所述方程(1)和(2)如下：

$\begin{array}{c}{E}^{inc}\left(r\right)=E\left(r\right)-[k_b^2+▿▿·]\underset{D_{inv}}{\int }g(r,r^{\prime }){\chi }^E\left(r^{\prime }\right)E\left(r^{\prime }\right){\mathrm{dr}}^{\prime }\\ +{\mathrm{j\omega \mu }}_b▿\times \underset{D_{inv}}{\int }g(r,r^{\prime }){\chi }^H\left(r^{\prime }\right)H\left(r^{\prime }\right){\mathrm{dr}}^{\prime }\end{array}---\left(1\right)$

$\begin{array}{c}{H}^{inc}\left(r\right)=H\left(r\right)-[k_b^2+▿▿·]\underset{D_{inv}}{\int }g(r,r^{\prime }){\chi }^H\left(r^{\prime }\right)H\left(r^{\prime }\right){\mathrm{dr}}^{\prime }\\ +{\mathrm{j\omega ϵ}}_b^{*}▿\times \underset{D_{inv}}{\int }g(r,r^{\prime }){\chi }^E\left(r^{\prime }\right)E\left(r^{\prime }\right){\mathrm{dr}}^{\prime }\end{array}---\left(2\right)$

其中，D_inv代表反演区域，r和r’为反演区域中的位置，k_b为背景介质的波数，ω为角频率，j为虚数符号；ε_b^*，μ_b分别是背景介质的复介电常数和磁导率；E^inc(r)和E(r)分别为入射电场和总电场，H^inc(r)和H(r)分别为入射磁场和总磁场；g(r,r’)是标量格林函，χ^E(r) 是电特性对比度，χ^H(r)是磁化率，定义为以下方程：

$g(r,r^{\prime })=\frac{e^{-{\mathrm{jk}}_b\left|r-r^{\prime }\right|}}{4\pi \left|r-r^{\prime }\right|}---\left(3\right)$

${\chi }^E\left(r\right)=\frac{{ϵ}^{*}\left(r\right)-{ϵ}_b^{*}}{{ϵ}_b^{*}}---\left(4\right)$

${\chi }^H\left(r\right)=\frac{\mu \left(r\right)-{\mu }_b}{{\mu }_b}---\left(5\right)$

其中，方程(4)和(5)中的ε^*(r)和μ(r)分别是人体复介电常数和磁导率，而复介电常数又可以用方程(6)表示：

$ϵ*\left(r\right)=ϵ\left(r\right)-j\frac{\sigma \left(r\right)}{\omega }---\left(6\right)$

其中，ε(r)，σ(r)分别是人体的介电常数，电导率；

②将第n-1次迭代的总电场E(r)和总磁场H(r)，代入方程(7)和(8)，进行反演计算，解出第n次迭代的χ^E；(n)(r)和χ^H；(n)(r)；

所述方程(7)和(8)如下：

$\begin{array}{c}{H}^{sca}\left(r\right)={\mathrm{j\omega ϵ}}_b^{*}▿\times \underset{D_{inv}}{\int }g(r,r^{\prime }){\chi }^E\left(r^{\prime }\right)E\left(r^{\prime }\right){\mathrm{dr}}^{\prime }\\ +[k_b^2+▿▿·]\underset{D_{inv}}{\int }g(r,r^{\prime }){\chi }^H\left(r^{\prime }\right)H\left(r^{\prime }\right){\mathrm{dr}}^{\prime }\end{array}---\left(7\right)$

$H_1^{+;sca}\left(r\right)=\frac{B_1^{+}\left(r\right)}{\mu \left(r\right)}-\frac{B_1^{+;inc}\left(r\right)}{{\mu }_b}---\left(8\right)$

其中，H^sca(r)散射磁场；H₁^+；sca(r)为正旋散射磁场，定义为H₁^+；sca(r)＝[H_x^sca(r)+jH_y^sca(r)]/2；H_x^sca(r)和H_y^sca(r)分别为散射磁场H^sca(r)的x分量与y分量；B₁⁺(r)为磁共振正旋场的磁通密度，B₁^+；inc(r)为背景介质下磁共振正旋场的磁通密度；

③将步骤②解出的χ^E；(n)(r)和χ^H；(n)(r)，代入方程(1)和(2)，进行正演计算，解出第n次迭代的E⁽ⁿ⁾(r)和H⁽ⁿ⁾(r)；

④重复步骤②和③，直到利用χ^E；(n)(r)，χ^H；(n)(r)代入方程(7)和(8)算出的B₁⁺(r)与测量值的相对残差在10％以内，则算法收敛，结束迭代。

所述迭代的算法可采用玻恩迭代算法(BIM)、变分玻恩迭代算法(VBIM)或变形玻恩迭代算法(DBIM)等。

在步骤2)第①和③步骤中，所述正演计算为己知电特性对比度χ^E(r)和磁化率χ^H(r)，利用方程(1)和(2)，结合稳定双共轭梯度算法(BCGS)和快速傅里叶变换(FFT)，计算出总电场E(r)和总磁场H(r)。

在步骤2)第②步骤中，所述反演计算为己知总电场E(r)和总磁场H(r)，利用方程(7)和(8)，结合共轭梯度算法(CG)和快速傅里叶变换(FFT)，计算出电特性对比度χ^E(r)和磁化率χ^H(r)。

3)利用方程(4),(5),(6)结合迭代结果中的电特性对比度χ^E(r)和磁化率χ^H(r),计算出ε(r)，σ(r)，μ(r)(或磁化率χ^H(r))的分布，输出结果，计算完毕。

本发明具体涉及电和磁体积分方程，结合玻恩迭代算法(BIM)或变分玻恩迭代算法(VBIM)或变形玻恩迭代算法(DBIM)、快速傅里叶变换(FFT)、稳定双共轭梯度算法(BCGS)以及共轭梯度算法(CG)，进行反复的正演和反演迭代至结果收敛，求解人体电磁特性参数(包括电导率σ，介电常数ε和磁导率μ)，实现人体电磁特性参数磁共振成像的方法，得到人体组织电磁特性参数分布的二维或三维图像，呈现组织内各区域电磁特性参数差异的同时，也可对其进行定量研究，可用于医学上的疾病研究和指导临床诊断、治疗。

本发明能有效抑制边界出现的错误、减小噪声影响、提高成像分辨率，结合MR信号的幅值信息和相位信息，利用快速傅里叶变换(FFT)加快重建速度、减少重建时间，结合玻恩迭代算法(BIM)或变分玻恩迭代算法(VBIM)或变形玻恩迭代算法(DBIM)，稳定双共轭梯度算法(BCGS)，共轭梯度算法(CG),同时反演电特性参数(电导率σ，介电常数ε)和磁特性参数(磁导率μ)。

本发明的有益技术效果如下：

提出一种基于快速体积分方程算法和磁共振的人体电磁特性反演方法，对于电特性成像，能够有效抑制边界出现的错误、减小噪声影响、提高成像分辨率，利用快速傅里叶变换(FFT)加快重建速度、减少重建时间，结合变分玻恩迭代算法(BIM)或变分玻恩迭代算法(VBIM)或变形玻恩迭代算法(DBIM)，稳定双共轭梯度算法(BCGS)，共轭梯度算法(CG)，同时反演电特性参数(电导率σ、介电常数ε)和磁特性参数(磁导率μ)。

附图说明

图1是反演出的人脑相对介电常数三维分布。

图2是反演出的人脑电导率三维分布。

图3是反演出的人脑磁化率三维分布。

图4是反演出的人脑相对介电常数分布X-Y截面。

图5是反演出的人脑电导率分布X-Y截面。

图6是反演出的人脑磁化率分布X-Y截面。

具体实施方式

本发明利用快速体积分方程算法和磁共振数据进行人体的二维或三维电磁特性反演，其中快速体积分方程算法中结合了玻恩迭代算法(BIM)或变分玻恩迭代算法(VBIM)或变形玻恩迭代算法(DBIM)，快速傅里叶变换(FFT)，稳定双共轭梯度算法(BCGS)，共轭梯度算法(CG)。

本实施例以变分玻恩迭代算法(VBIM)为例子进行解释。具体实施方式如下：

1)测量实验数据。利用磁共振的影像技术测量人体，得到B₁⁺场。B₁⁺场是磁共振射频场的正旋场。

2)进行迭代计算，包括以下步骤：

①假定初始的电特性对比度χ^E；(0)(r)和磁化率χ^H；(0)(r),进行正演计算，解出E⁽⁰⁾(r)和H⁽⁰⁾(r)。

②利用第n-1次迭代的总电场E(r)和总磁场H(r)，进行反演计算，解出第n次迭代的χ^E；(n)(r)和χ^H；(n)(r)。

③再由步骤②解出的χ^E；(n)(r)和χ^H；(n)(r)，进行正演计算，解出第n次迭代的E⁽ⁿ⁾(r)和H⁽ⁿ⁾(r)。

④重复步骤②和③，直到χ^E；(n)(r)，χ^H；(n)(r)代入方程(14)和(15)算出的B₁⁺(r)与测量值的相对残差在10％以内，则算法收敛，结束迭代。

步骤①和③所述正演计算为己知电特性对比度χ^E(r)和磁化率χ^H(r)，利用以下方程计算总电场E(r)和总磁场H(r):

$E^{inc}\left(r\right)=E\left(r\right)+j\omega [1+\frac{1}{k_b^2}▿▿·]A\left(r\right)+\frac{▿\times F\left(r\right)}{{ϵ}_b^{*}}---\left(1\right)$

$H^{inc}\left(r\right)=H\left(r\right)+j\omega [1+\frac{1}{k_b^2}▿▿·]F\left(r\right)-\frac{▿\times A\left(r\right)}{{\mu }_b}---\left(2\right)$

$A\left(r\right)={\mathrm{j\omega \mu }}_b{ϵ}_b^{*}\underset{D_{inv}}{\int }g(r,r^{\prime }){\chi }^E\left(r^{\prime }\right)E\left(r^{\prime }\right){\mathrm{dr}}^{\prime }---\left(3\right)$

$F\left(r\right)={\mathrm{j\omega \mu }}_b{ϵ}_b^{*}\underset{D_{inv}}{\int }g(r,r^{\prime }){\chi }^H\left(r^{\prime }\right)H\left(r^{\prime }\right){\mathrm{dr}}^{\prime }---\left(4\right)$

其中，D_inv代表反演区域，r和r’为反演区域中的位置，k_b为背景介质的波数，ω为角频率，j为虚数符号。ε_b^*，μ_b分别是背景介质的复介电常数和磁导率。E^inc(r)和E(r)分别为入射电场和总电场，H^inc(r)和H(r)分别为入射磁场和总磁场。A(r)和F(r)分别为磁势矢量和电势矢量。g(r,r’)是标量格林函数，χ^E(r)是电特性对比度，χ^H(r)是磁化率，定义为以下方程：

$g(r,r^{\prime })=\frac{e^{-{\mathrm{jk}}_b\left|r-r^{\prime }\right|}}{4\pi \left|r-r^{\prime }\right|}---\left(5\right)$

${\chi }^E\left(r\right)=\frac{{ϵ}^{*}\left(r\right)-{ϵ}_b^{*}}{{ϵ}_b^{*}}---\left(6\right)$

${\chi }^H\left(r\right)=\frac{\mu \left(r\right)-{\mu }_b}{{\mu }_b}---\left(7\right)$

其中，方程(4)和(5)中的ε^*(r)和μ(r)分别是人体复介电常数和磁导率，而复介电常数又可以用方程(6)来表示：

${ϵ}^{*}\left(r\right)=ϵ\left(r\right)-j\frac{\sigma \left(r\right)}{\omega }---\left(8\right)$

其中，ε(r)，σ(r)分别是人体的介电常数，电导率。

便于计算，利用冲激函数作为基函数和测试函数，则方程(1)，(2)，(3)和(4)被离散化为下列方程：

$E_{j,k,l}^{inc}=E_{j,k,l}+j\omega [1+\frac{1}{k_b^2}▿▿·]A_{j,k,l}+\frac{▿\times F_{j,k,l}}{{ϵ}_b^{*}}---\left(9\right)$

$H_{j,k,l}^{inc}=H_{j,k,l}+j\omega [1+\frac{1}{k_b^2}▿▿·]F_{j,k,l}-\frac{▿\times A_{j,k,l}}{{\mu }_b}---\left(10\right)$

$A_{j,k,l}={\mathrm{j\omega \mu }}_b{ϵ}_b^{*}\Delta V\underset{j^{\prime }=1}{\overset{J}{\Sigma }}\underset{k^{\prime }=1}{\overset{K}{\Sigma }}\underset{l^{\prime }=1}{\overset{L}{\Sigma }}{g}_{j-j^{\prime },k-k^{\prime },l-l^{\prime }}{\chi }_{j^{\prime },k^{\prime },l^{\prime }}^E{E}_{j^{\prime },k^{\prime },l^{\prime }}---\left(11\right)$

$F_{j,k,l}={\mathrm{j\omega \mu }}_b{ϵ}_b^{*}\Delta V\underset{j^{\prime }=1}{\overset{J}{\Sigma }}\underset{k^{\prime }=1}{\overset{K}{\Sigma }}\underset{l^{\prime }=1}{\overset{L}{\Sigma }}{g}_{j-j^{\prime },k-k^{\prime },l-l^{\prime }}{\chi }_{j^{\prime },k^{\prime },l^{\prime }}^H{E}_{j^{\prime },k^{\prime },l^{\prime }}---\left(12\right)$

其中，离散的单元格中心为：

$r_{j,k,l}=(x_j,y_k,z_l)=\left(\right(j-\frac{1}{2})\Delta x,(k-\frac{1}{2})\Delta y,(1-\frac{1}{2}\left)\Delta z\right)$

j∈[1,J],k∈[1,K],l∈[1,L]

Δx、Δy和Δz分别为单元格在x、y和z方向的大小。j，k，l分别为离散后x方向，y方向和z方向的位置。J，K，L为x方向，y方向和z方向的离散点数。

而方程(11)和(12)中的ΔV为单元格体积，定义为ΔV＝ΔxΔyΔz。

由于，标量格林函数g(r,r’)具有平移不变性，所以方程(11)和(12)可以利用快速傅里叶变换(FFT)进行快速计算。此方法可以节省大量时间以及内存，使计算更加快速。

离散化后，可以将方程写成以下形式：

其中，

$W=\left(\begin{array}{c}E\\ H\end{array}\right),W^{inc}=\left(\begin{array}{c}{E}^{inc}\\ H^{inc}\end{array}\right)$

L代表的是线性操作。

直接解方程(13)会花大量时间和内存，所以在本发明当中利用稳定双共轭梯度算法(BCGS)和快速傅里叶变换(FFT)计算方程(13)得出总电场E(r)和总磁场H(r)。

步骤②所述反演计算为己知总电场E(r)和总磁场H(r)，利用以下方程计算电特性对比度χ^E(r)和磁化率χ^H(r):

$H^{sca}\left(r\right)={\mathrm{j\omega ϵ}}_b^{*}\underset{D_{inv}}{\int }{G}_b^H(r,r^{\prime }){\chi }^E\left(r^{\prime }\right){\mathrm{dr}}^{\prime }{\mathrm{dr}}^{\prime }+k_b^2\underset{D_{inv}}{\int }{G}_b^E(r,r^{\prime }){\chi }^H\left(r^{\prime }\right)H\left(r^{\prime }\right){\mathrm{dr}}^{\prime }---\left(14\right)$

$H_1^{+;sca}\left(r\right)=\frac{B_1^{+}\left(r\right)}{\mu \left(r\right)}-\frac{B_1^{+;inc}\left(r\right)}{{\mu }_b}---\left(15\right)$

其中，H^sca(r)散射磁场。H₁^+；sca(r)为正旋散射磁场，定义为H₁^+；sca(r)＝[H_x^sca(r)+jH_y^sca(r)]/2，其中，H_x^sca(r)和H_y^sca(r)分别为散射磁场H^sca(r)的x分量与y分量。B₁⁺(r)为磁共振正旋场的磁通密度，B₁^+；inc(r)为背景介质下磁共振正旋场的磁通密度。G_b^E(r,r’)和G_b^H(r,r’)分别为电并矢格林函数和磁并格林函数。

本发明中利用变分玻恩迭代来计算物体的电磁特性分布。所以方程(14)写成以下形式：

${\mathrm{\delta H}}^{sca}\left(r\right)={\mathrm{j\omega ϵ}}_b^{*}\underset{D_{inv}}{\int }{G}_b^H(r,r^{\prime }){\mathrm{\delta \chi }}^E\left(r^{\prime }\right)E\left(r^{\prime }\right){\mathrm{dr}}^{\prime }+k_b^2\underset{D_{inv}}{\int }{G}_b^E(r,r^{\prime }){\mathrm{\delta \chi }}^H\left(r^{\prime }\right)H\left(r^{\prime }\right){\mathrm{dr}}^{\prime }---\left(16\right)$

其中，

δH^sca(r)＝H^sca(r)-H^sca；(n-1)(r)>

代表的是测量值与n-1次迭代的散射磁场的差值。

δχ^E(r)＝χ^E；(n)(r)-χ^E；(n-1)(r)>

δχ^H(r)＝χ^H；(n)(r)-χ^H；(n-1)(r)>

代表的是第n次迭代与第n-1次迭代计算的电特性对比度χ^E(r)和磁化率χ^H(r)的差值。

为了便于计算，利用冲激函数作为基函数和测试函数，则方程(16)写为以下形式：

$\begin{array}{c}{\mathrm{\delta H}}_{j,k,l}^{sca}={\mathrm{j\omega ϵ}}_b^{*}\underset{j^{\prime }=1}{\overset{J}{\Sigma }}\underset{k^{\prime }=1}{\overset{K}{\Sigma }}\underset{l^{\prime }=1}{\overset{L}{\Sigma }}{G}_{j-j^{\prime },k-k^{\prime },l-l}^{\prime }{{\mathrm{\delta \chi }}^E}_{j^{\prime },k^{\prime },l^{\prime }}{E}_{j^{\prime },k^{\prime },l^{\prime }}\\ +k_b^2\underset{j^{\prime }=1}{\overset{J}{\Sigma }}\underset{k^{\prime }=1}{\overset{K}{\Sigma }}\underset{l^{\prime }=1}{\overset{L}{\Sigma }}{G}_{j-j^{\prime },k-k^{\prime },l-l}^E{{\mathrm{\delta \chi }}^H}_{j^{\prime },k^{\prime },l^{\prime }}{H}_{j^{\prime },k^{\prime },l^{\prime }}\end{array}---\left(20\right)$

其中，电并矢格林函数G_b^E(r,r’)和磁并格林函数G_b^H(r,r’)都具有平移不变性，所以也可以利用快速傅里叶变换(FFT)来进行快速计算。

当利用方程(20)计算δχ^E(r)和δχ^H(r)时，如果直接用矩阵方式进行计算的话，花费的时间与内存过于庞大，所以本发明利用共轭梯度算法(CG)，结合快速傅里叶变换来计算δχ^E(r)和δχ^H(r)，然后利用方程(18)和(19)算出χ^E；(n)(r)和χ^H；(n)(r)。

3)利用方程(6),(7),(8)结合迭代结果中的电特性对比度χ^E(r)和磁化率χ^H(r),计算出ε(r)，σ(r)，μ(r)(或磁化率χ^H(r))的分布，输出结果，计算完毕。

表1

组织相对介电常数电导率(S/m)磁化率脑干980.5110小脑1170.7190灰质980.5116×10^-8白质680.2913×10^-8头骨170.0600下丘脑980.5112×10^-8眼睛851.000舌头750.6520

人脑设置参数如表1,频率为63.87MHz。