基于瞬时线化的电特性断层成像方法、系统、设备及介质

1.本发明属于核磁共振成像技术领域，尤其涉及一种基于瞬时线化的电特性断层成像方法、系统、设备及介质。


背景技术：

2.对于高场核磁共振成像系统而言，频率会随着场强的变大而变大，波长会随着频率的变大而变小，直至与扫描物体的波长在同一水平。这种状态会引发共振，会使被测物体与射频场相互作用变强，可以利用这种现象来实现成像。上述基于高场磁共振的方法称为介电特性断层成像(mr-ert)，是一种新的核磁共振量化成像方法。
3.生物组织中存在大量水和少量无机盐等电特性较强的物质，而不同组织间物质的分布又各不相同。该特性决定了生物组织是一种电特性较强的物体，可以利用此特性来成像。不同于传统断层成像方法，介电特性断层成像是一种定量测量，这意味着可以量化地进行分析，对癌症的早起分析等一类传统特性成像较难解决的问题，介电特性断层成像有着先天优势。除此之外，介电特性会影响脑内电信号的传播，并与脑组织的类型、活动、功能、衰老、病理改变等信息高度，相关因此mr-ept有望为研究脑认知与脑疾病提供全新的研究工具。人体内介电特性分布还是计算特定吸收率(sar)的基础，可以用来评估微波场和射频场下的人体安全。热效应是超高场磁共振成像系统一直需要面对的一个很严重的安全问题，而对于人体组织的介电特性的准确成像能够有助于sar的研究，推动超高场磁共振系统能够更安全地应用。
4.介电特性断层成像分为两个步骤，第一步是射频场信息的采集。第二步是以麦克斯韦方程为基础，得到介电特性断层成像的方程，然后构造求解电特性分布的数学模型，并根据已采集的射频场信息作为模型的输入，来重建电特性的分布。其中第二步是电特性成像算法所需要研究的重点，如何构建正确的数学模型来更加准确地重建电特性分布。
5.目前，用于电特性重建的主流算法，大多数是忽略一些次要问题，或者对方程本身的复杂参数进行简化。有一种经典算法的前提条件是电特性本身分布是均匀的，这种方法被称为helmholtz ept。这种方法的缺陷非常明显，就是实际的电特性分布不可能全是均匀的，这意味着方法本身的前提就是不符合现实条件的。这会使得均匀假设失效，从而导致严重的计算误差，会在非均匀区产生很多错误数据。而且，亥姆霍兹方程自身是一个高阶非线性非齐次的偏微分复数域方程，这意味着如果不进行简化，方程本身很难进行求解。


技术实现要素：

6.本发明的目的在于提供一种基于瞬时线化的电特性断层成像方法，能够在尽量少简化或者更改问题模型的前提下，更加准确地得到电特性断层成像。
7.本发明是通过以下技术方案实现的：
8.一种基于瞬时线化的电特性断层成像方法，包括以下步骤：
9.s1、测量磁共振射频场，基于射频场计算得到发射场；
10.s2、利用数值方法计算得到发射场的点数据信息，点数据信息包括发射场在空间中每一点处沿x轴、y轴以及z轴方向的一阶导数和二阶导数；
11.s3、基于发射场和发射场的点数据信息，根据麦克斯韦方程组，计算得到空间中每一点处的电容率初值和电导率初值，并由空间中每一点处的电容率初值和电导率初值确定约束项k、r1与r2；
12.s4、基于发射场、发射场的点数据信息以及约束项k、r1与r2，根据麦克斯韦方程组结合牛顿迭代法，得到迭代方程，迭代方程如公式(1)所示：
[0013][0014]
式中，un表示第n代的u，f(u)为由发射场b的点数据信息确定的函数，α为阻尼因子向量；
[0015]
s5、基于空间中每一处的电容率初值和电导率初值，确定迭代方程的初值，将初值代入迭代方程中进行迭代计算，经过若干次迭代后得到空间中每一点处的电容率和电导率；
[0016]
s6、基于空间中每一点处的电容率和电导率输出电特性断层成像图。
[0017]
进一步地，在迭代方程中，f(u)如公式(2)所示：
[0018]
f(u)＝au+cu2ꢀꢀ
(2)
[0019]
式中，u2为一列向量，u2的每一项的通式为a为矩阵，如公式(3)所示：
[0020][0021]
式中，变量如公式(4)所示，变量γ如公式(5)所示：
[0022][0023][0024]
c为矩阵，如公式(6)所示：
[0025][0026]
进一步地，在迭代方程中，阻尼因子向量α的计算方式如下：
[0027]
构建一个nx*ny*nz的矩阵q，q每一点的值如公式(7)所示：
[0028][0029]
式中，σm和εm均为参数；
[0030]
设一矩阵v，对于v中每一点v,均有在q中对应的以v坐标为中心的m*m矩阵t，v的值如公式(23)所示：
[0031]
v＝min(t)
ꢀꢀ
(8)
[0032]
将矩阵v一一映射到阻尼因子向量α中即得到阻尼因子向量α。
[0033]
进一步地，由空间中每一点处的电容率初值和电导率初值确定约束项k、r1与r2的步骤包括：
[0034]
对于空间中每一点处的电容率初值和电导率初值，判断电容率初值和电导率初值沿各方向导数是否为0；
[0035]
若是，则将电容率初值和电导率初值记为目标电容率初值和目标电导率初值，目标电容率初值和目标电导率初值不参与迭代，并根据目标电容率初值和目标电导率初值构造约束项k；
[0036]
若否，则将电容率初值和电导率初值记为迭代电容率初值和迭代电导率初值；
[0037]
约束项r1为规模为1*(n
x
*ny*nz)的列向量，约束项r1中每一项的通式如公式(9)所示：
[0038][0039]
式中，λ
r1
为约束项r1的系数，ω为核磁共振设备对应的拉莫尔进动频率，n0与(x0,y0,z0)满足如公式(10)所示的关系：
[0040]
n0＝nz*ny*x0+nz*y0+z0ꢀꢀ
(10)
[0041]fr1
如式(11)所示：
[0042][0043]
式中，θr为参数；
[0044]
常数项θ
t0
如公式(12)所示：
[0045][0046]
式中，常数p为目标电容率初值的数量；
[0047]
约束项r2为一规模为1*(n
x
*ny*nz)的列向量，约束项r2如公式(13)所示：
[0048][0049]
式中，λ
r2
为约束项r2的系数，为拉普拉斯算子，u为一规模为1*(n
x
*ny*nz)的列向量，u的每一项通式如公式(14)所示：
[0050][0051]
进一步地，基于空间中每一处的电容率初值和电导率初值，确定迭代方程的初值，将初值代入迭代方程中进行迭代计算，经过若干次迭代后得到空间中每一点处的电容率和电导率的步骤包括：
[0052]
s51、对迭代电容率初值和迭代电导率初值进行预处理，基于预处理后的迭代电容率初值和迭代电导率初值，根据公式(15)计算迭代方程的初值u0：
[0053]
u0＝σ0+iωε0ꢀꢀ
(15)
[0054]
式中，σ0为迭代电导率初，ε0为迭代电容率初值；
[0055]
s52、将初值u0代入迭代方程中，得到新的初值，判断新的初值中的各项平均值是否小于预设的迭代阈值；
[0056]
s53、若否，则重复步骤s52；
[0057]
s54、若是，则根据迭代后新的初值、目标电容率初值和目标电导率初值得到空间中每一点处的电容率和电导率。
[0058]
进一步地，利用数值方法计算得到发射场的点数据信息的步骤包括：
[0059]
根据公式(16)计算发射场在空间中每一点处沿x轴方向的一阶导数
[0060][0061]
根据公式(17)计算发射场在空间中每一点处沿y轴方向的一阶导数
[0062][0063]
根据公式(18)计算发射场在空间中每一点处沿z轴方向的一阶导数
[0064][0065]
根据公式(19)计算发射场在空间中每一点处沿x轴方向的二阶导数
[0066][0067]
根据公式(20)计算发射场在空间中每一点处沿y轴方向的二阶导数
[0068][0069]
根据公式(21)计算发射场在空间中每一点处沿z轴方向的二阶导数
[0070][0071]
其中，dx、dy、dz分别为空间离散化后的单位网格在x轴、y轴以及z轴方向上的长度。
[0072]
进一步地，测量磁共振射频场，基于射频场计算得到发射场的步骤包括：
[0073]
测量磁共振射频场沿空间坐标系x轴、y轴以及z轴方向的分量b
x
、by、bz；
[0074]
根据公式(22)计算发射场b：
[0075]
[0076]
式中，i为复数单位。
[0077]
本发明还提供了一种基于瞬时线化的电特性断层成像系统，包括：
[0078]
测量模块，用于测量磁共振射频场，基于射频场计算得到发射场；
[0079]
第一计算模块，用于利用数值方法计算得到发射场的点数据信息，点数据信息包括发射场在空间中每一点处沿x轴、y轴以及z轴方向的一阶导数和二阶导数；
[0080]
第二计算模块，用于基于发射场和发射场的点数据信息，根据麦克斯韦方程组，计算得到空间中每一点处的电容率初值和电导率初值，并由空间中每一点处的电容率初值和电导率初值确定约束项k、r1与r2；
[0081]
得到模块，用于基于发射场、发射场的点数据信息以及约束项k、r1与r2，根据麦克斯韦方程组结合牛顿迭代法，得到迭代方程，迭代方程如公式(1)所示：
[0082][0083]
式中，un表示第n代的u，f(u)为由发射场b的点数据信息确定的函数，α为阻尼因子向量；
[0084]
迭代模块，用于基于空间中每一处的电容率初值和电导率初值，确定迭代方程的初值，将初值代入迭代方程中进行迭代计算，经过若干次迭代后得到空间中每一点处的电容率和电导率；
[0085]
输出模块，用于基于空间中每一点处的电容率和电导率输出电特性断层成像图。
[0086]
本发明还公开了一种电子设备，包括存储器和处理器，存储器存储有计算机程序，处理器执行计算机程序时实现上述任一项方法的步骤。
[0087]
本发明还公开了一种计算机可读存储介质，其上存储有计算机程序，计算机程序被处理器执行时实现上述任一项的方法的步骤。
[0088]
相比于现有技术，本发明的有益效果为：将原始麦克斯韦方程和亥姆霍兹方程推导得到的高阶非线性非齐次偏微分方程通过改良的牛顿迭代法实现瞬时线化得到迭代方程，并根据物理模型和现实情况加入约束项，构成完备的具有初值的迭代方程，得到一个精确度更高，更符合真实状况的电特性分布解，从而能够在尽量少简化或者更改问题模型的前提下，更加准确地得到电特性断层成像。
附图说明
[0089]
图1为本发明基于瞬时线化的电特性断层成像方法的步骤流程图；
[0090]
图2为本发明基于瞬时线化的电特性断层成像方法与helmholtz ept方法、改良crept方法在duke人脑肿瘤仿真模型上的重建电导率(上)与电容率(下)图像的结果对比；
[0091]
图3为本发明基于瞬时线化的电特性断层成像方法与helmholtz ept方法，改良crept方法在另一例duke人脑肿瘤仿真模型上的重建电导率(上)与电容率(下)图像的结果；
[0092]
图4为本发明基于瞬时线化的电特性断层成像系统的模块示意图；
[0093]
图5为本发明电子设备一实施例的结构示意图；
[0094]
图6为本发明计算机可读存储介质一实施例的结构示意框图。
具体实施方式
[0095]
为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。通常在此处附图中描述和示出的本发明实施例的组件可以以各种不同的配置来布置和设计。
[0096]
因此，以下对在附图中提供的本发明的实施例的详细描述并非旨在限制要求保护的本发明的范围，而是仅仅表示本发明的选定实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
[0097]
应注意到：相似的标号和字母在下面的附图中表示类似项，因此，一旦某一项在一个附图中被定义，则在随后的附图中不需要对其进行进一步定义和解释。同时，在本发明的描述中，术语“第一”、“第二”等仅用于区分描述，而不能理解为指示或暗示相对重要性。
[0098]
需要说明的是，在本文中，诸如第一和第二等之类的关系术语仅仅用来将一个实体或者操作与另一个实体或操作区分开来，而不一定要求或者暗示这些实体或操作之间存在任何这种实际的关系或者顺序。而且，术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法、物品或者设备不仅包括那些要素，而且还包括没有明确列出的其他要素，或者是还包括为这种过程、方法、物品或者设备所固有的要素。在没有更多限制的情况下，由语句“包括一个
……”
限定的要素，并不排除在包括所述要素的过程、方法、物品或者设备中还存在另外的相同要素。
[0099]
在本发明的描述中，需要说明的是，术语“上”、“下”、“内”、“外”等指示的方位或位置关系为基于附图所示的方位或位置关系，或者是该发明产品使用时惯常摆放的方位或位置关系，仅是为了便于描述本发明和简化描述，而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作，因此不能理解为对本发明的限制。
[0100]
请参阅图1，图1为本发明基于瞬时线化的电特性断层成像方法的步骤流程图。一种基于瞬时线化的电特性断层成像方法，包括以下步骤：
[0101]
s1、测量磁共振射频场，基于射频场计算得到发射场；
[0102]
s2、利用数值方法计算得到发射场的点数据信息，点数据信息包括发射场在空间中每一点处沿x轴、y轴以及z轴方向的一阶导数和二阶导数；
[0103]
s3、基于发射场和发射场的点数据信息，根据麦克斯韦方程组，计算得到空间中每一点处的电容率初值和电导率初值，并由空间中每一点处的电容率初值和电导率初值确定约束项k、r1与r2；
[0104]
s4、基于发射场、发射场的点数据信息以及约束项k、r1与r2，根据麦克斯韦方程组结合牛顿迭代法，得到迭代方程，迭代方程如公式(1)所示：
[0105][0106]
式中，un表示第n代的u，即第n次迭代后得到的迭代值，f(u)为由发射场b的点数据信息确定的函数，α为阻尼因子向量；
[0107]
s5、基于空间中每一处的电容率初值和电导率初值，确定迭代方程的初值，将初值代入迭代方程中进行迭代计算，经过若干次迭代后得到空间中每一点处的电容率和电导
率；
[0108]
s6、基于空间中每一点处的电容率和电导率输出电特性断层成像图。
[0109]
在上述步骤s1中，用核磁共振设备配合正交线圈以及特定b1扫描序列可以得到正交线圈测量空间的射频场b，即得到磁共振射频场b。然后基于得到的射频场b，计算出发射场b。
[0110]
进一步地，在步骤s1中，测量磁共振射频场，基于射频场计算得到发射场的步骤包括：
[0111]
s11、测量磁共振射频场沿空间坐标系x轴、y轴以及z轴方向的分量b
x
、by、bz；
[0112]
s12、根据公式(22)计算发射场b：
[0113][0114]
式中，i为复数单位。
[0115]
在上述步骤s11至步骤s12中，由仪器测得的射频场b唯一矢量，其值为(b
x
,by,bz)，其中，b
x
、by、bz分别为射频场b沿空间坐标系x轴、y轴以及z轴方向的分量。而发射场b由公式计算。
[0116]
在上述步骤s2中，数值方法可以采用有限差分方法，基于发射场b，通过有限差分方法计算得到发射场在空间中每一点处沿x轴、y轴以及z轴方向的一阶导数和二阶导数。
[0117]
进一步地，利用数值方法计算得到发射场的点数据信息的步骤包括：
[0118]
根据公式(16)计算发射场在空间中每一点处沿x轴方向的一阶导数
[0119][0120]
根据公式(17)计算发射场在空间中每一点处沿y轴方向的一阶导数
[0121][0122]
根据公式(18)计算发射场在空间中每一点处沿z轴方向的一阶导数
[0123][0124]
根据公式(19)计算发射场在空间中每一点处沿x轴方向的二阶导数
[0125][0126]
根据公式(20)计算发射场在空间中每一点处沿y轴方向的二阶导数
[0127][0128]
根据公式(21)计算发射场在空间中每一点处沿z轴方向的二阶导数
[0129][0130]
其中，dx、dy、dz分别为空间离散化后的单位网格在x轴、y轴以及z轴方向上的长度，i为角标，表示空间中特定一点。
[0131]
在上述步骤s3中，基于发射场和发射场的点数据信息，根据mrept(magnetic resonance electrical properties，人体组织电特性磁共振断层成像)的原理，计算得到空间中每一点处的电容率初值ε0和电导率初值σ0，然后根据预设的处理规则对计算得到空间中每一点处的电容率初值ε0和电导率初值σ0进行处理，确定约束项k、r1与r2。
[0132]
进一步地，在步骤s3中，由空间中每一点处的电容率初值和电导率初值确定约束项k、r1与r2的步骤包括：
[0133]
s31、对于空间中每一点处的电容率初值和电导率初值，判断电容率初值和电导率初值沿各方向导数是否为0；
[0134]
s32、若是，则将电容率初值和电导率初值记为目标电容率初值和目标电导率初值，目标电容率初值和目标电导率初值不参与迭代，并根据目标电容率初值和目标电导率初值构造约束项k；
[0135]
s33、若否，则将电容率初值和电导率初值记为迭代电容率初值和迭代电导率初值；
[0136]
s34、约束项r1为规模为1*(n
x
*ny*nz)的列向量，约束项r1中每一项的通式如公式(9)所示：
[0137][0138]
式中，λ
r1
为约束项r1的系数，ω为核磁共振设备对应的拉莫尔进动频率，n0与(x0,y0,z0)满足如公式(10)所示的关系：
[0139]
n0＝nz*ny*x0+nz*y0+z0ꢀꢀ
(10)
[0140]fr1
如式(11)所示：
[0141][0142]
式中，θr为参数，一般情况下取值0.2；
[0143]
常数项θ
t0
如公式(12)所示：
[0144][0145]
式中，常数p为目标电容率初值的数量；
[0146]
s35、约束项r2为一规模为1*(n
x
*ny*nz)的列向量，约束项r2如公式(13)所示：
[0147][0148]
式中，λ
r2
为约束项r2的系数，为拉普拉斯算子，u为一规模为1*(n
x
*ny*nz)的列向量，u的每一项通式如公式(14)所示：
[0149][0150]
在上述步骤s31至步骤s33中，参考helmholtz ept方法，假设电特性参数在空间中分布均匀，即空间中每一点处的电特性参数沿各方向导数均为0，这样可以将偏微分方程转化为线性方程，便于得到一个用于迭代的初值，其中电特性参数包括电容率和电导率。而在
真实空间分布中确实有一部分电容率初值ε0和电导率初值σ0是满足沿各方向的导数等于0的假设条件的，该部分电容率初值ε0和电导率初值σ0已经非常准确，没有必要再次进行求解，可以从待求量变为已知量来约束迭代方程本身，这样可以降低后续迭代计算的计算复杂度，进而降低计算量。因此，基于helmholtz ept方法的假设前提，对于空间中每一点处的电容率初值ε0和电导率初值σ0，分析电容率初值ε0和电导率初值σ0沿各方向导数是否为0，若是，则说明该点处的电容率初值ε0和电导率初值σ0是正确的，将该点处的电容率初值ε0和电导率初值σ0记为目标电容率初值ε
t0
和目标电导率初值σ
t0
，不用参与迭代，可以从待求量变为已知量k来约束迭代方程，约束项k为向量，由目标电容率初值ε0和电导率初值σ0构成，具体的构成方式由偏微分方程和矩阵方程的方程变形分解原理计算得到，在此不再赘述；若否，则说明该点处的电容率初值ε0和电导率初值σ0是不合理的，需要参与迭代，将该点处的电容率初值ε0和电导率初值σ0记为迭代电容率初值ε
d0
和电导率初值σ
d0
，作为待求量代入迭代方程中进行求解。参与迭代的迭代电容率初值ε
d0
和电导率初值σ
d0
的规模为n
x
*ny*nz，n
x
、ny、nz表示x，y，z方向上像素的个数，那么整个空间的像素个数即为n
x
*ny*nz。
[0151]
在上述步骤s34中，目标电导率初值的数量与目标电容率初值的数量相同，因此常数p也可以为目标电导率初值的数量，例如，在步骤s31和步骤s32中筛选出1000个合理的目标电容率初值和目标电导率初值，即目标电容率初值和目标电导率初值的数量均为1000，因此常数p为1000。
[0152]
在上述步骤s4中，基于发射场和发射场的点数据信息，将原始麦克斯韦方程和亥姆霍兹方程推导得到的高阶非线性非齐次偏微分方程通过改良的牛顿迭代法实现瞬时线化得到迭代公式，并根据物理模型和现实中电容率初值ε0和电导率初值σ0的情况加入约束项k、r1与r2，得到迭代方程
[0153]
进一步地，在迭代方程中，f(u)如公式(2)所示：
[0154]
f(u)＝au+cu2ꢀꢀ
(2)
[0155]
式中，u2为一列向量，u2的每一项的通式为a为由发射场和发射场的点数据信息确定的矩阵，如公式(3)所示：
[0156][0157]
式中，变量如公式(4)所示，变量γ如公式(5)所示：
[0158][0159][0160]
c为由发射场和发射场的点数据信息确定的矩阵，如公式(6)所示：
[0161][0162]
进一步地，在迭代方程中，阻尼因子向量α的计算方式如下：
[0163]
构建一个n
x
*ny*nz的矩阵q，q每一点的值如公式(7)所示：
[0164][0165]
式中，σm和εm均为参数，(xi,yi,zi)与角标i一一对应，i表示第i个点；
[0166]
设一矩阵v，对于v中每一点v,均有在q中对应的以v坐标为中心的m*m矩阵t，v的值如公式(23)所示：
[0167]
v＝min(t)
ꢀꢀ
(8)
[0168]
将矩阵v一一映射到阻尼因子向量α中即得到阻尼因子向量α。
[0169]
矩阵v的元素v对应空间中的每一个点，矩阵t相当于一个扫描窗口，其行数和列数m可以根据实际情况确定，优选地，m为5。而对于5*5矩阵t，每一个元素v，在矩阵q中以对应的v坐标的点为中心点，以中心点和中心点各方向上的两个点构建一个5*5矩阵t，在构建的5*5矩阵t中找到最小的值作为元素v的值，以此不断得到每一个元素v的值，进而得到矩阵v。并且由上述过程可知，矩阵v在每次迭代时会改变，因此每次迭代后需要重新计算矩阵v，以重新计算阻尼因子向量α。
[0170]
在上述步骤s5中，在迭代法中有一个合理的初值以加快收敛速度，因此通过根据空间中每一处的电容率初值和电导率初值，确定迭代方程的初值可以提高收敛速度。而本发明中磁共振使用的设备是正交线圈，故发射场b在z方向上分量可忽略不计，故所涉及原始方程如公式(23)所示：
[0171][0172]
式中，i表示虚数，ω为核磁共振设备对应的拉莫尔进动频率，μ0为磁导率，b
′
x
为yz平面的发射场b沿x轴的偏导数，b
′y为xz平面的发射场b沿y轴的偏导数，b
′z为xy平面的发射场b沿z轴的偏导数；
[0173]
根据helmholtz ept方法，假设电特性参数在空间中分布均匀，即空间中每一点处的电特性参数沿各方向导数均为0，使得公式(23)可以简化为如公式(24)所示：
[0174][0175]
由公式(24)推导可得到公式(25)：
[0176]
[0177]
进而得到迭代方程所需要的初值u0＝σ0+iωε0，构造的初值u0时，是从一个由空间中每一处的电容率初值和电导率初值构成的矩阵变换为一个向量u0，具有一一对应的关系，将初值u0代入迭代方程中进行迭代计算，经过若干次迭代后得到空间中每一点处的电容率和电导率。
[0178]
进一步地，基于空间中每一处的电容率初值和电导率初值，确定迭代方程的初值，将初值代入迭代方程中进行迭代计算，经过若干次迭代后得到空间中每一点处的电容率和电导率的步骤包括：
[0179]
s51、对迭代电容率初值和迭代电导率初值进行预处理，基于预处理后的迭代电容率初值和迭代电导率初值，根据公式(15)计算迭代方程的初值u0：
[0180]
u0＝σ0+iωε0ꢀꢀ
(15)
[0181]
式中，σ0为迭代电导率初，ε0为迭代电容率初值；
[0182]
s52、将初值u0代入迭代方程中，得到新的初值，判断新的初值中的各项平均值是否小于预设的迭代阈值；
[0183]
s53、若否，则重复步骤s52；
[0184]
s54、若是，则根据迭代后新的初值、目标电容率初值和目标电导率初值得到空间中每一点处的电容率和电导率。
[0185]
在上述步骤s51至步骤s54中，由步骤s31至步骤s33筛选出需要进行迭代的迭代电容率初值ε
d0
和电导率初值σ
d0
，根据公式(15)构成初值u0，将初值u0代入迭代方程中，得到新的初值u1，新的初值u1中的各项平均值是否小于预设的迭代阈值，判断f(u1)的均值是否小迭代阈值，若否，则根据新的初值u1代入迭代方程，得到新的初值u2，以此不断进行迭代，直到某一代的初值ui中的各项平均值小于预设的迭代阈值，则结束迭代，得到ui。由于构件初值u0时，是从一个矩阵变换为一个向量，具有一一对应的关系，因此可以事先记录每个迭代电容率初值ε
d0
和电导率初值σ
d0
的映射信息，及空间中对应的坐标，同样事先记录由步骤s31至步骤s33筛选出来不需要参与迭代的目标电容率初值ε
t0
和目标电导率初值σ
t0
的映射信息，在计算得到ui后，基于ui、目标电容率初值ε
t0
和目标电导率初值σ
t0
，根据事先记录的坐标信息还原成矩阵形式，即可得到空间中每一点处的电容率值ε和电导率值σ。
[0186]
在上述步骤s6中，空间中每一点处的电容率值ε和电导率值σ即为图像上每一点的信息，即可输出电特性断层成像图。
[0187]
以下是本发明基于瞬时线化的电特性断层成像方法(以下简称本发明)与helmholtz ept和改良的crept方法的结果对比：
[0188]
参考图为duke人脑肿瘤仿真模型，评估指标选择峰值信噪比以及结构相似性来比较在本例中本发明与helmholtz ept和改良的crept方法的结果。峰值信噪比与结构相似性的值越高，说明重建出的电特性断层成像图与参考图像越接近，表1是本发明与helmholtz ept和改良的crept方法得出的结果的量化指标。
[0189]
表1
[0190]
方法名称峰值信噪比结构相似性helmholtz eptu8.61510.6825改良creptu15.52950.8178本发明u18.24230.8684
[0191]
从表1可见，在本例中，无论是峰值信噪比还是结构相似性本发明均优于helmholtz ept和改良的crept方法。
[0192]
另外结合请结合参阅图2和图3，图2为本发明基于瞬时线化的电特性断层成像方法与helmholtz ept方法、改良crept方法在duke人脑肿瘤仿真模型上的重建电导率(上)与电容率(下)图像的结果对比，图3为本发明基于瞬时线化的电特性断层成像方法与helmholtz ept方法，改良crept方法在另一例duke人脑肿瘤仿真模型上的重建电导率(上)与电容率(下)图像的结果。由图2和图3可知，本发明图视觉上更接近于参考图。因此，本发明得到的电特性断层成像图的效果要显著好于helmholtz ept和改良的crept方法。
[0193]
请结合参阅图4，图4为本发明基于瞬时线化的电特性断层成像系统的模块示意图。本发明还提供了一种基于瞬时线化的电特性断层成像系统，包括：
[0194]
测量模块1，用于测量磁共振射频场，基于射频场计算得到发射场；
[0195]
第一计算模块2，用于利用数值方法计算得到发射场的点数据信息，点数据信息包括发射场在空间中每一点处沿x轴、y轴以及z轴方向的一阶导数和二阶导数；
[0196]
第二计算模块3，用于基于发射场和发射场的点数据信息，根据麦克斯韦方程组，计算得到空间中每一点处的电容率初值和电导率初值，并由空间中每一点处的电容率初值和电导率初值确定约束项k、r1与r2；
[0197]
得到模块4，用于基于发射场、发射场的点数据信息以及约束项k、r1与r2，根据麦克斯韦方程组结合牛顿迭代法，得到迭代方程，迭代方程如公式(1)所示：
[0198][0199]
式中，un表示第n代的u，即第n次迭代后得到的迭代值，f(u)为由发射场b的点数据信息确定的函数，α为阻尼因子向量；
[0200]
迭代模块5，用于基于空间中每一处的电容率初值和电导率初值，确定迭代方程的初值，将初值代入迭代方程中进行迭代计算，经过若干次迭代后得到空间中每一点处的电容率和电导率；
[0201]
输出模块6，用于基于空间中每一点处的电容率和电导率输出电特性断层成像图。
[0202]
在上述测量模块1中，用核磁共振设备配合正交线圈以及特定b1扫描序列可以得到正交线圈测量空间的射频场b，即得到磁共振射频场b。然后基于得到的射频场b，计算出发射场b。
[0203]
进一步地，测量模块1包括：
[0204]
测量子模块，测量磁共振射频场沿空间坐标系x轴、y轴以及z轴方向的分量b
x
、by、bz；
[0205]
第一计算子模块，用于根据公式(22)计算发射场b：
[0206][0207]
式中，i为复数单位。
[0208]
在上述测量子模块至第一计算子模块中，由仪器测得的射频场b唯一矢量，其值为(b
x
,by,bz)，其中，b
x
、by、bz分别为射频场b沿空间坐标系x轴、y轴以及z轴方向的分量。而发射场b由公式计算。
[0209]
数值方法可以采用有限差分方法，基于发射场b，第一计算模块2通过有限差分方法计算得到发射场在空间中每一点处沿x轴、y轴以及z轴方向的一阶导数和二阶导数。
[0210]
进一步地，第一计算模块2包括：
[0211]
第二计算子模块，用于根据公式(16)计算发射场在空间中每一点处沿x轴方向的一阶导数
[0212][0213]
第三计算子模块，用于根据公式(17)计算发射场在空间中每一点处沿y轴方向的一阶导数
[0214][0215]
第四计算子模块，用于根据公式(18)计算发射场在空间中每一点处沿z轴方向的一阶导数
[0216][0217]
第五计算子模块，用于根据公式(19)计算发射场在空间中每一点处沿x轴方向的二阶导数
[0218][0219]
第六计算子模块，用于根据公式(20)计算发射场在空间中每一点处沿y轴方向的二阶导数
[0220][0221]
第七计算子模块，用于根据公式(21)计算发射场在空间中每一点处沿z轴方向的二阶导数
[0222][0223]
其中，dx、dy、dz分别为空间离散化后的单位网格在x轴、y轴以及z轴方向上的长度，i为角标，表示空间中特定一点。
[0224]
第二计算模块3基于发射场和发射场的点数据信息，根据mrept(magnetic resonance electrical properties，人体组织电特性磁共振断层成像)的原理，计算得到空间中每一点处的电容率初值ε0和电导率初值σ0，然后根据预设的处理规则对计算得到空间中每一点处的电容率初值ε0和电导率初值σ0进行处理，确定约束项k、r1与r2。
[0225]
进一步地，第二计算模块3包括：
[0226]
判断子模块，用于对于空间中每一点处的电容率初值和电导率初值，判断电容率初值和电导率初值沿各方向导数是否为0；
[0227]
第一标记子模块，用于若判断子模块的判断为是，则将电容率初值和电导率初值记为目标电容率初值和目标电导率初值，目标电容率初值和目标电导率初值不参与迭代，并根据目标电容率初值和目标电导率初值构造约束项k；
[0228]
第二标记子模块，用于若判断子模块的判断为否，则将电容率初值和电导率初值记为迭代电容率初值和迭代电导率初值；
[0229]
第一约束项模块，用于约束项r1为规模为1*(n
x
*ny*nz)的列向量，约束项r1中每一项的通式如公式(9)所示：
[0230][0231]
式中，λ
r1
为约束项r1的系数，ω为核磁共振设备对应的拉莫尔进动频率，n0与(x0,y0,z0)满足如公式(10)所示的关系：
[0232]
n0＝nz*ny*x0+nz*y0+z0ꢀꢀ
(10)
[0233]fr1
如式(11)所示：
[0234][0235]
式中，θr为参数，一般情况下取值0.2；
[0236]
常数项θ
t0
如公式(12)所示：
[0237][0238]
式中，常数p为目标电容率初值的数量；
[0239]
第一约束项模块，用于约束项r2为一规模为1*(n
x
*ny*nz)的列向量，约束项r2如公式(13)所示：
[0240][0241]
式中，λ
r2
为约束项r2的系数，为拉普拉斯算子，u为一规模为1*(n
x
*ny*nz)的列向量，u的每一项通式如公式(14)所示：
[0242][0243]
参考helmholtz ept方法，假设电特性参数在空间中分布均匀，即空间中每一点处的电特性参数沿各方向导数均为0，这样可以将偏微分方程转化为线性方程，便于得到一个用于迭代的初值，其中电特性参数包括电容率和电导率。而在真实空间分布中确实有一部分电容率初值ε0和电导率初值σ0是满足沿各方向的导数等于0的假设条件的，该部分电容率初值ε0和电导率初值σ0已经非常准确，没有必要再次进行求解，可以从待求量变为已知量来约束迭代方程本身，这样可以降低后续迭代计算的计算复杂度，进而降低计算量。因此，基于helmholtz ept方法的假设前提，判断子模块对于空间中每一点处的电容率初值ε0和电导率初值σ0，分析电容率初值ε0和电导率初值σ0沿各方向导数是否为0，若是，则说明该点处的电容率初值ε0和电导率初值σ0是正确的，第一标记模块将该点处的电容率初值ε0和电导率初值σ0记为目标电容率初值ε
t0
和目标电导率初值σ
t0
，不用参与迭代，可以从待求量变为已知量k来约束迭代方程，约束项k为向量，由目标电容率初值ε0和电导率初值σ0构成，具体的构成方式由偏微分方程和矩阵方程的方程变形分解原理计算得到，在此不再赘述；若否，则说明该点处的电容率初值ε0和电导率初值σ0是不合理的，需要参与迭代，第二标记模块将
该点处的电容率初值ε0和电导率初值σ0记为迭代电容率初值ε
d0
和电导率初值σ
d0
，作为待求量代入迭代方程中进行求解。参与迭代的迭代电容率初值ε
d0
和电导率初值σ
d0
的规模为n
x
*ny*nz，n
x
、ny、nz表示x，y，z方向上像素的个数，那么整个空间的像素个数即为n
x
*ny*nz。
[0244]
在上述第一约束项模块中，目标电导率初值的数量与目标电容率初值的数量相同，因此常数p也可以为目标电导率初值的数量，例如，由判断子模块和第一标记子模块筛选出1000个合理的目标电容率初值和目标电导率初值，即目标电容率初值和目标电导率初值的数量均为1000，因此常数p为1000。
[0245]
得到模块4基于发射场和发射场的点数据信息，将原始麦克斯韦方程和亥姆霍兹方程推导得到的高阶非线性非齐次偏微分方程通过改良的牛顿迭代法实现瞬时线化得到迭代公式，并根据物理模型和现实中电容率初值ε0和电导率初值σ0的情况加入约束项k、r1与r2，得到迭代方程
[0246]
进一步地，在迭代方程中，f(u)如公式(2)所示：
[0247]
f(u)＝au+cu2ꢀꢀ
(2)
[0248]
式中，u2为一列向量，u2的每一项的通式为a为由发射场和发射场的点数据信息确定的矩阵，如公式(3)所示：
[0249][0250]
式中，变量如公式(4)所示，变量γ如公式(5)所示：
[0251][0252][0253]
c为由发射场和发射场的点数据信息确定的矩阵，如公式(6)所示：
[0254][0255]
进一步地，在迭代方程中，阻尼因子向量α的计算方式如下：
[0256]
构建一个n
x
*ny*nz的矩阵q，q每一点的值如公式(7)所示：
[0257][0258]
式中，σm和εm均为参数，(xi,yi,zi)与角标i一一对应，i表示第i个点；
[0259]
设一矩阵v，对于v中每一点v,均有在q中对应的以v坐标为中心的m*m矩阵t，v的值如公式(23)所示：
[0260]
v＝min(t)
ꢀꢀ
(8)
[0261]
将矩阵v一一映射到阻尼因子向量α中即得到阻尼因子向量α。
[0262]
矩阵v的元素v对应空间中的每一个点，矩阵t相当于一个扫描窗口，其行数和列数m可以根据实际情况确定，优选地，m为5。而对于5*5矩阵t，每一个元素v，在矩阵q中以对应的v坐标的点为中心点，以中心点和中心点各方向上的两个点构建一个5*5矩阵t，在构建的5*5矩阵t中找到最小的值作为元素v的值，以此不断得到每一个元素v的值，进而得到矩阵v。并且由上述过程可知，矩阵v在每次迭代时会改变，因此每次迭代后需要重新计算矩阵v，以重新计算阻尼因子向量α。
[0263]
在上述迭代模块5中，在迭代法中有一个合理的初值以加快收敛速度，因此通过根据空间中每一处的电容率初值和电导率初值，确定迭代方程的初值可以提高收敛速度。而本发明中磁共振使用的设备是正交线圈，故发射场b在z方向上分量可忽略不计，故所涉及原始方程如公式(23)所示：
[0264][0265]
式中，i表示虚数，ω为核磁共振设备对应的拉莫尔进动频率，μ0为磁导率，b
′
x
为yz平面的发射场b沿x轴的偏导数，b
′y为xz平面的发射场b沿y轴的偏导数，b
′z为xy平面的发射场b沿z轴的偏导数；
[0266]
根据helmholtz ept方法，假设电特性参数在空间中分布均匀，即空间中每一点处的电特性参数沿各方向导数均为0，使得公式(23)可以简化为如公式(24)所示：
[0267][0268]
由公式(24)推导可得到公式(25)：
[0269][0270]
进而得到迭代方程所需要的初值u0＝σ0+iωε0，构造的初值u0时，是从一个由空间中每一处的电容率初值和电导率初值构成的矩阵变换为一个向量u0，具有一一对应的关系，将初值u0代入迭代方程中进行迭代计算，经过若干次迭代后得到空间中每一点处的电容率和电导率。
[0271]
进一步地，迭代模块5包括：
[0272]
处理子模块，用于对迭代电容率初值和迭代电导率初值进行预处理，基于预处理后的迭代电容率初值和迭代电导率初值，根据公式(15)计算迭代方程的初值u0：
[0273]
u0＝σ0+iωε0ꢀꢀ
(15)
[0274]
式中，σ0为迭代电导率初，ε0为迭代电容率初值；
[0275]
代入模块，用于将初值u0代入迭代方程中，得到新的初值，判断新的初值中的各项平均值是否小于预设的迭代阈值；
[0276]
重复子模块，用于若否，则重复步骤s52；
[0277]
得到子模块，用于若是，则根据迭代后新的初值、目标电容率初值和目标电导率初值得到空间中每一点处的电容率和电导率。
[0278]
在上述处理子模块至得到子模块中，由判断子模块、第一标记子模块和第二标记子模块筛选出需要进行迭代的迭代电容率初值ε
d0
和电导率初值σ
d0
，根据公式(15)构成初值u0，将初值u0代入迭代方程中，得到新的初值u1，新的初值u1中的各项平均值是否小于预设的迭代阈值，判断f(u1)的均值是否小迭代阈值，若否，则根据新的初值u1代入迭代方程，得到新的初值u2，以此不断进行迭代，直到某一代的初值ui中的各项平均值小于预设的迭代阈值，则结束迭代，得到ui。由于构件初值u0时，是从一个矩阵变换为一个向量，具有一一对应的关系，因此可以事先记录每个迭代电容率初值ε
d0
和电导率初值σ
d0
的映射信息，及空间中对应的坐标，同样事先记录由步骤s31至步骤s33筛选出来不需要参与迭代的目标电容率初值ε
t0
和目标电导率初值σ
t0
的映射信息，在计算得到ui后，基于ui、目标电容率初值ε
t0
和目标电导率初值σ
t0
，根据事先记录的坐标信息还原成矩阵形式，即可得到空间中每一点处的电容率值ε和电导率值σ。
[0279]
在上述输出模块6中，空间中每一点处的电容率值ε和电导率值σ即为图像上每一点的信息，即可输出电特性断层成像图。
[0280]
请结合参阅图5，图5为本发明电子设备一实施例的结构示意框图。本发明一实施例还提出一种电子设备1001，包括存储器1003和处理器1002，存储器1003存储有计算机程序1004，处理器1002执行计算机程序1004时实现上述任一项基于瞬时线化的电特性断层成像方法的步骤，包括：s1、测量磁共振射频场，基于射频场计算得到发射场；s2、利用数值方法计算得到发射场的点数据信息，点数据信息包括发射场在空间中每一点处沿x轴、y轴以及z轴方向的一阶导数和二阶导数；s3、基于发射场和发射场的点数据信息，根据麦克斯韦方程组，计算得到空间中每一点处的电容率初值和电导率初值，并由空间中每一点处的电容率初值和电导率初值确定约束项k、r1与r2；s4、基于发射场、发射场的点数据信息以及约束项k、r1与r2，根据麦克斯韦方程组结合牛顿迭代法，得到迭代方程，根据麦克斯韦方程组结合牛顿迭代法，得到迭代方程s5、基于空间中每一处的电容率初值和电导率初值，确定迭代方程的初值，将初值代入迭代方程中进行迭代计算，经过若干次迭代后得到空间中每一点处的电容率和电导率；s6、基于空间中每一点处的电容率和电导率输出电特性断层成像图。
[0281]
请结合参阅图6，图6为本发明计算机可读存储介质一实施例的结构示意框图。本发明一实施例还提供一种计算机可读存储介质2001，其上存储有计算机程序1004，计算机程序1004被处理器1002执行时实现上述任一项基于瞬时线化的电特性断层成像方法的步骤，包括：s1、测量磁共振射频场，基于射频场计算得到发射场；s2、利用数值方法计算得到发射场的点数据信息，点数据信息包括发射场在空间中每一点处沿x轴、y轴以及z轴方向的一阶导数和二阶导数；s3、基于发射场和发射场的点数据信息，根据麦克斯韦方程组，计算得到空间中每一点处的电容率初值和电导率初值，并由空间中每一点处的电容率初值和电
导率初值确定约束项k、r1与r2；s4、基于发射场、发射场的点数据信息以及约束项k、r1与r2，根据麦克斯韦方程组结合牛顿迭代法，得到迭代方程根据麦克斯韦方程组结合牛顿迭代法，得到迭代方程s5、基于空间中每一处的电容率初值和电导率初值，确定迭代方程的初值，将初值代入迭代方程中进行迭代计算，经过若干次迭代后得到空间中每一点处的电容率和电导率；s6、基于空间中每一点处的电容率和电导率输出电特性断层成像图。
[0282]
相比于现有技术，本发明的有益效果为：将原始麦克斯韦方程和亥姆霍兹方程推导得到的高阶非线性非齐次偏微分方程通过改良的牛顿迭代法实现瞬时线化得到迭代方程，并根据物理模型和现实情况加入约束项，构成完备的具有初值的迭代方程，得到一个精确度更高，更符合真实状况的电特性分布解，从而能够在尽量少简化或者更改问题模型的前提下，更加准确地得到电特性断层成像。
[0283]
本领域普通技术人员可以理解实现上述实施例方法中的全部或部分流程，是可以通过计算机程序来指令相关的硬件来完成，所述的计算机程序可存储于一非易失性计算机可读取存储介质中，该计算机程序在执行时，可包括如上述各方法的实施例的流程。其中，本技术所提供的和实施例中所使用的对存储器、存储、数据库或其它介质的任何引用，均可包括非易失性和/或易失性存储器。非易失性存储器可以包括只读存储器(rom)、可编程rom(prom)、电可编程rom(eprom)、电可擦除可编程rom(eeprom)或闪存。易失性存储器可包括随机存取存储器(ram)或者外部高速缓冲存储器。作为说明而非局限，ram以多种形式可得，诸如静态ram(sram)、动态ram(dram)、同步dram(sdram)、双速据率sdram(ssrsdram)、增强型sdram(esdram)、同步链路(synchlink)dram(sldram)、存储器总线(rambus)直接ram(rdram)、直接存储器总线动态ram(drdram)、以及存储器总线动态ram(rdram)等。
[0284]
需要说明的是，在本文中，术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性的包含，从而使得包括一系列要素的过程、装置、物品或者方法不仅包括那些要素，而且还包括没有明确列出的其他要素，或者是还包括为这种过程、装置、物品或者方法所固有的要素。在没有更多限制的情况下，由语句“包括一个
……”
限定的要素，并不排除在包括该要素的过程、装置、物品或者方法中还存在另外的相同要素。
[0285]
本发明并不局限于上述实施方式，如果对本发明的各种改动或变形不脱离本发明的精神和范围，倘若这些改动和变形属于本发明的权利要求和等同技术范围之内，则本发明也意图包含这些改动和变形。