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

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

Info

Publication number
CN115808650A
CN115808650A CN202211343568.4A CN202211343568A CN115808650A CN 115808650 A CN115808650 A CN 115808650A CN 202211343568 A CN202211343568 A CN 202211343568A CN 115808650 A CN115808650 A CN 115808650A
Authority
CN
China
Prior art keywords
equation
permittivity
conductivity
initial value
point
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Pending
Application number
CN202211343568.4A
Other languages
English (en)
Inventor
陈武凡
王兆年
阮国辉
夏灵
路利军
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Southern Medical University
Original Assignee
Southern Medical University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Southern Medical University filed Critical Southern Medical University
Priority to CN202211343568.4A priority Critical patent/CN115808650A/zh
Publication of CN115808650A publication Critical patent/CN115808650A/zh
Pending legal-status Critical Current

Links

Images

Landscapes

  • Magnetic Resonance Imaging Apparatus (AREA)

Abstract

本发明提供了一种基于瞬时线化的电特性断层成像方法、系统、设备及介质，其中方法包括：测量磁共振射频场，基于射频场计算得到发射场；利用数值方法计算得到发射场的点数据信息；基于发射场和发射场的点数据信息，根据麦克斯韦方程组，得到空间中每一点处的电容率初值和电导率初值，并由电容率初值和电导率初值确定约束项K、R1与R2；基于发射场、发射场的点数据信息以及约束项K、R1与R2，根据麦克斯韦方程组结合牛顿迭代法得到迭代方程；基于电容率初值和电导率初值确定迭代方程的初值，将初值代入迭代方程中进行迭代计算，经过若干次迭代后得到空间中每一点处的电容率和电导率；基于空间中每一点处的电容率和电导率输出电特性断层成像图。

Description

基于瞬时线化的电特性断层成像方法、系统、设备及介质
技术领域
本发明属于核磁共振成像技术领域，尤其涉及一种基于瞬时线化的电特性断层成像方法、系统、设备及介质。
背景技术
对于高场核磁共振成像系统而言，频率会随着场强的变大而变大，波长会随着频率的变大而变小，直至与扫描物体的波长在同一水平。这种状态会引发共振，会使被测物体与射频场相互作用变强，可以利用这种现象来实现成像。上述基于高场磁共振的方法称为介电特性断层成像(MR-ERT)，是一种新的核磁共振量化成像方法。
生物组织中存在大量水和少量无机盐等电特性较强的物质，而不同组织间物质的分布又各不相同。该特性决定了生物组织是一种电特性较强的物体，可以利用此特性来成像。不同于传统断层成像方法，介电特性断层成像是一种定量测量，这意味着可以量化地进行分析，对癌症的早起分析等一类传统特性成像较难解决的问题，介电特性断层成像有着先天优势。除此之外，介电特性会影响脑内电信号的传播，并与脑组织的类型、活动、功能、衰老、病理改变等信息高度，相关因此MR-EPT有望为研究脑认知与脑疾病提供全新的研究工具。人体内介电特性分布还是计算特定吸收率(SAR)的基础，可以用来评估微波场和射频场下的人体安全。热效应是超高场磁共振成像系统一直需要面对的一个很严重的安全问题，而对于人体组织的介电特性的准确成像能够有助于SAR的研究，推动超高场磁共振系统能够更安全地应用。
介电特性断层成像分为两个步骤，第一步是射频场信息的采集。第二步是以麦克斯韦方程为基础，得到介电特性断层成像的方程，然后构造求解电特性分布的数学模型，并根据已采集的射频场信息作为模型的输入，来重建电特性的分布。其中第二步是电特性成像算法所需要研究的重点，如何构建正确的数学模型来更加准确地重建电特性分布。
目前，用于电特性重建的主流算法，大多数是忽略一些次要问题，或者对方程本身的复杂参数进行简化。有一种经典算法的前提条件是电特性本身分布是均匀的，这种方法被称为Helmholtz EPT。这种方法的缺陷非常明显，就是实际的电特性分布不可能全是均匀的，这意味着方法本身的前提就是不符合现实条件的。这会使得均匀假设失效，从而导致严重的计算误差，会在非均匀区产生很多错误数据。而且，亥姆霍兹方程自身是一个高阶非线性非齐次的偏微分复数域方程，这意味着如果不进行简化，方程本身很难进行求解。
发明内容
本发明的目的在于提供一种基于瞬时线化的电特性断层成像方法，能够在尽量少简化或者更改问题模型的前提下，更加准确地得到电特性断层成像。
本发明是通过以下技术方案实现的：
一种基于瞬时线化的电特性断层成像方法，包括以下步骤：
S1、测量磁共振射频场，基于射频场计算得到发射场；
S2、利用数值方法计算得到发射场的点数据信息，点数据信息包括发射场在空间中每一点处沿x轴、y轴以及z轴方向的一阶导数和二阶导数；
S3、基于发射场和发射场的点数据信息，根据麦克斯韦方程组，计算得到空间中每一点处的电容率初值和电导率初值，并由空间中每一点处的电容率初值和电导率初值确定约束项K、R1与R2；
S4、基于发射场、发射场的点数据信息以及约束项K、R1与R2，根据麦克斯韦方程组结合牛顿迭代法，得到迭代方程，迭代方程如公式(1)所示：
Figure BDA0003917416640000021
式中，un表示第n代的u，f(u)为由发射场B的点数据信息确定的函数，α为阻尼因子向量；
S5、基于空间中每一处的电容率初值和电导率初值，确定迭代方程的初值，将初值代入迭代方程中进行迭代计算，经过若干次迭代后得到空间中每一点处的电容率和电导率；
S6、基于空间中每一点处的电容率和电导率输出电特性断层成像图。
进一步地，在迭代方程中，f(u)如公式(2)所示：
f(u)＝Au+Cu2 (2)
式中，u2为一列向量，u2的每一项的通式为
Figure BDA0003917416640000022
A为矩阵，如公式(3)所示：
Figure BDA0003917416640000023
式中，变量
Figure BDA0003917416640000024
如公式(4)所示，变量γ如公式(5)所示：
Figure BDA0003917416640000031
Figure BDA0003917416640000032
C为矩阵，如公式(6)所示：
Figure BDA0003917416640000033
进一步地，在迭代方程中，阻尼因子向量α的计算方式如下：
构建一个nx*ny*nz的矩阵Q，Q每一点的值如公式(7)所示：
Figure BDA0003917416640000034
式中，σm和εm均为参数；
设一矩阵V，对于V中每一点v,均有在Q中对应的以v坐标为中心的m*m矩阵T，v的值如公式(23)所示：
v＝min(T) (8)
将矩阵V一一映射到阻尼因子向量α中即得到阻尼因子向量α。
进一步地，由空间中每一点处的电容率初值和电导率初值确定约束项K、R1与R2的步骤包括：
对于空间中每一点处的电容率初值和电导率初值，判断电容率初值和电导率初值沿各方向导数是否为0；
若是，则将电容率初值和电导率初值记为目标电容率初值和目标电导率初值，目标电容率初值和目标电导率初值不参与迭代，并根据目标电容率初值和目标电导率初值构造约束项K；
若否，则将电容率初值和电导率初值记为迭代电容率初值和迭代电导率初值；
约束项R1为规模为1*(nx*ny*nz)的列向量，约束项R1中每一项的通式如公式(9)所示：
Figure BDA0003917416640000035
式中，λR1为约束项R1的系数，
Figure BDA0003917416640000041
ω为核磁共振设备对应的拉莫尔进动频率，n0与(x0,y0,z0)满足如公式(10)所示的关系：
n0＝nz*ny*x0+nz*y0+z0 (10)
fR1如式(11)所示：
Figure BDA0003917416640000042
式中，θr为参数；
常数项θt0如公式(12)所示：
Figure BDA0003917416640000043
式中，常数p为目标电容率初值的数量；
约束项R2为一规模为1*(nx*ny*nz)的列向量，约束项R2如公式(13)所示：
Figure BDA0003917416640000044
式中，λR2为约束项R2的系数，
Figure BDA0003917416640000047
为拉普拉斯算子，u为一规模为1*(nx*ny*nz)的列向量，u的每一项通式如公式(14)所示：
Figure BDA0003917416640000045
进一步地，基于空间中每一处的电容率初值和电导率初值，确定迭代方程的初值，将初值代入迭代方程中进行迭代计算，经过若干次迭代后得到空间中每一点处的电容率和电导率的步骤包括：
S51、对迭代电容率初值和迭代电导率初值进行预处理，基于预处理后的迭代电容率初值和迭代电导率初值，根据公式(15)计算迭代方程的初值u0：
u0＝σ0+iωε0 (15)
式中，σ0为迭代电导率初，ε0为迭代电容率初值；
S52、将初值u0代入迭代方程中，得到新的初值，判断新的初值中的各项平均值是否小于预设的迭代阈值；
S53、若否，则重复步骤S52；
S54、若是，则根据迭代后新的初值、目标电容率初值和目标电导率初值得到空间中每一点处的电容率和电导率。
进一步地，利用数值方法计算得到发射场的点数据信息的步骤包括：
根据公式(16)计算发射场在空间中每一点处沿x轴方向的一阶导数
Figure BDA0003917416640000046
Figure BDA0003917416640000051
根据公式(17)计算发射场在空间中每一点处沿y轴方向的一阶导数
Figure BDA0003917416640000052
Figure BDA0003917416640000053
根据公式(18)计算发射场在空间中每一点处沿z轴方向的一阶导数
Figure BDA0003917416640000054
Figure BDA0003917416640000055
根据公式(19)计算发射场在空间中每一点处沿x轴方向的二阶导数
Figure BDA0003917416640000056
Figure BDA0003917416640000057
根据公式(20)计算发射场在空间中每一点处沿y轴方向的二阶导数
Figure BDA0003917416640000058
Figure BDA0003917416640000059
根据公式(21)计算发射场在空间中每一点处沿z轴方向的二阶导数
Figure BDA00039174166400000510
Figure BDA00039174166400000511
其中，dx、dy、dz分别为空间离散化后的单位网格在x轴、y轴以及z轴方向上的长度。
进一步地，测量磁共振射频场，基于射频场计算得到发射场的步骤包括：
测量磁共振射频场沿空间坐标系x轴、y轴以及z轴方向的分量Bx、By、Bz；
根据公式(22)计算发射场B：
Figure BDA00039174166400000512
式中，i为复数单位。
本发明还提供了一种基于瞬时线化的电特性断层成像系统，包括：
测量模块，用于测量磁共振射频场，基于射频场计算得到发射场；
第一计算模块，用于利用数值方法计算得到发射场的点数据信息，点数据信息包括发射场在空间中每一点处沿x轴、y轴以及z轴方向的一阶导数和二阶导数；
第二计算模块，用于基于发射场和发射场的点数据信息，根据麦克斯韦方程组，计算得到空间中每一点处的电容率初值和电导率初值，并由空间中每一点处的电容率初值和电导率初值确定约束项K、R1与R2；
得到模块，用于基于发射场、发射场的点数据信息以及约束项K、R1与R2，根据麦克斯韦方程组结合牛顿迭代法，得到迭代方程，迭代方程如公式(1)所示：
Figure BDA00039174166400000513
式中，un表示第n代的u，f(u)为由发射场B的点数据信息确定的函数，α为阻尼因子向量；
迭代模块，用于基于空间中每一处的电容率初值和电导率初值，确定迭代方程的初值，将初值代入迭代方程中进行迭代计算，经过若干次迭代后得到空间中每一点处的电容率和电导率；
输出模块，用于基于空间中每一点处的电容率和电导率输出电特性断层成像图。
本发明还公开了一种电子设备，包括存储器和处理器，存储器存储有计算机程序，处理器执行计算机程序时实现上述任一项方法的步骤。
本发明还公开了一种计算机可读存储介质，其上存储有计算机程序，计算机程序被处理器执行时实现上述任一项的方法的步骤。
相比于现有技术，本发明的有益效果为：将原始麦克斯韦方程和亥姆霍兹方程推导得到的高阶非线性非齐次偏微分方程通过改良的牛顿迭代法实现瞬时线化得到迭代方程，并根据物理模型和现实情况加入约束项，构成完备的具有初值的迭代方程，得到一个精确度更高，更符合真实状况的电特性分布解，从而能够在尽量少简化或者更改问题模型的前提下，更加准确地得到电特性断层成像。
附图说明
图1为本发明基于瞬时线化的电特性断层成像方法的步骤流程图；
图2为本发明基于瞬时线化的电特性断层成像方法与Helmholtz EPT方法、改良crEPT方法在Duke人脑肿瘤仿真模型上的重建电导率(上)与电容率(下)图像的结果对比；
图3为本发明基于瞬时线化的电特性断层成像方法与Helmholtz EPT方法，改良crEPT方法在另一例Duke人脑肿瘤仿真模型上的重建电导率(上)与电容率(下)图像的结果；
图4为本发明基于瞬时线化的电特性断层成像系统的模块示意图；
图5为本发明电子设备一实施例的结构示意图；
图6为本发明计算机可读存储介质一实施例的结构示意框图。
具体实施方式
为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。通常在此处附图中描述和示出的本发明实施例的组件可以以各种不同的配置来布置和设计。
因此，以下对在附图中提供的本发明的实施例的详细描述并非旨在限制要求保护的本发明的范围，而是仅仅表示本发明的选定实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
应注意到：相似的标号和字母在下面的附图中表示类似项，因此，一旦某一项在一个附图中被定义，则在随后的附图中不需要对其进行进一步定义和解释。同时，在本发明的描述中，术语“第一”、“第二”等仅用于区分描述，而不能理解为指示或暗示相对重要性。
需要说明的是，在本文中，诸如第一和第二等之类的关系术语仅仅用来将一个实体或者操作与另一个实体或操作区分开来，而不一定要求或者暗示这些实体或操作之间存在任何这种实际的关系或者顺序。而且，术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法、物品或者设备不仅包括那些要素，而且还包括没有明确列出的其他要素，或者是还包括为这种过程、方法、物品或者设备所固有的要素。在没有更多限制的情况下，由语句“包括一个……”限定的要素，并不排除在包括所述要素的过程、方法、物品或者设备中还存在另外的相同要素。
在本发明的描述中，需要说明的是，术语“上”、“下”、“内”、“外”等指示的方位或位置关系为基于附图所示的方位或位置关系，或者是该发明产品使用时惯常摆放的方位或位置关系，仅是为了便于描述本发明和简化描述，而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作，因此不能理解为对本发明的限制。
请参阅图1，图1为本发明基于瞬时线化的电特性断层成像方法的步骤流程图。一种基于瞬时线化的电特性断层成像方法，包括以下步骤：
S1、测量磁共振射频场，基于射频场计算得到发射场；
S2、利用数值方法计算得到发射场的点数据信息，点数据信息包括发射场在空间中每一点处沿x轴、y轴以及z轴方向的一阶导数和二阶导数；
S3、基于发射场和发射场的点数据信息，根据麦克斯韦方程组，计算得到空间中每一点处的电容率初值和电导率初值，并由空间中每一点处的电容率初值和电导率初值确定约束项K、R1与R2；
S4、基于发射场、发射场的点数据信息以及约束项K、R1与R2，根据麦克斯韦方程组结合牛顿迭代法，得到迭代方程，迭代方程如公式(1)所示：
Figure BDA0003917416640000071
式中，un表示第n代的u，即第n次迭代后得到的迭代值，f(u)为由发射场B的点数据信息确定的函数，α为阻尼因子向量；
S5、基于空间中每一处的电容率初值和电导率初值，确定迭代方程的初值，将初值代入迭代方程中进行迭代计算，经过若干次迭代后得到空间中每一点处的电容率和电导率；
S6、基于空间中每一点处的电容率和电导率输出电特性断层成像图。
在上述步骤S1中，用核磁共振设备配合正交线圈以及特定B1扫描序列可以得到正交线圈测量空间的射频场B，即得到磁共振射频场B。然后基于得到的射频场B，计算出发射场B。
进一步地，在步骤S1中，测量磁共振射频场，基于射频场计算得到发射场的步骤包括：
S11、测量磁共振射频场沿空间坐标系x轴、y轴以及z轴方向的分量Bx、By、Bz；
S12、根据公式(22)计算发射场B：
Figure BDA0003917416640000081
式中，i为复数单位。
在上述步骤S11至步骤S12中，由仪器测得的射频场B唯一矢量，其值为(Bx,By,Bz)，其中，Bx、By、Bz分别为射频场B沿空间坐标系x轴、y轴以及z轴方向的分量。而发射场B由公式
Figure BDA0003917416640000082
计算。
在上述步骤S2中，数值方法可以采用有限差分方法，基于发射场B，通过有限差分方法计算得到发射场在空间中每一点处沿x轴、y轴以及z轴方向的一阶导数和二阶导数。
进一步地，利用数值方法计算得到发射场的点数据信息的步骤包括：
根据公式(16)计算发射场在空间中每一点处沿x轴方向的一阶导数
Figure BDA0003917416640000083
Figure BDA0003917416640000084
根据公式(17)计算发射场在空间中每一点处沿y轴方向的一阶导数
Figure BDA0003917416640000085
Figure BDA0003917416640000086
根据公式(18)计算发射场在空间中每一点处沿z轴方向的一阶导数
Figure BDA0003917416640000087
Figure BDA0003917416640000088
根据公式(19)计算发射场在空间中每一点处沿x轴方向的二阶导数
Figure BDA0003917416640000089
Figure BDA00039174166400000810
根据公式(20)计算发射场在空间中每一点处沿y轴方向的二阶导数
Figure BDA00039174166400000811
Figure BDA00039174166400000812
根据公式(21)计算发射场在空间中每一点处沿z轴方向的二阶导数
Figure BDA00039174166400000813
Figure BDA00039174166400000814
其中，dx、dy、dz分别为空间离散化后的单位网格在x轴、y轴以及z轴方向上的长度，i为角标，表示空间中特定一点。
在上述步骤S3中，基于发射场和发射场的点数据信息，根据MREPT(magneticresonance electrical properties，人体组织电特性磁共振断层成像)的原理，计算得到空间中每一点处的电容率初值ε0和电导率初值σ0，然后根据预设的处理规则对计算得到空间中每一点处的电容率初值ε0和电导率初值σ0进行处理，确定约束项K、R1与R2。
进一步地，在步骤S3中，由空间中每一点处的电容率初值和电导率初值确定约束项K、R1与R2的步骤包括：
S31、对于空间中每一点处的电容率初值和电导率初值，判断电容率初值和电导率初值沿各方向导数是否为0；
S32、若是，则将电容率初值和电导率初值记为目标电容率初值和目标电导率初值，目标电容率初值和目标电导率初值不参与迭代，并根据目标电容率初值和目标电导率初值构造约束项K；
S33、若否，则将电容率初值和电导率初值记为迭代电容率初值和迭代电导率初值；
S34、约束项R1为规模为1*(nx*ny*nz)的列向量，约束项R1中每一项的通式如公式(9)所示：
Figure BDA0003917416640000091
式中，λR1为约束项R1的系数，
Figure BDA0003917416640000092
ω为核磁共振设备对应的拉莫尔进动频率，n0与(x0,y0,z0)满足如公式(10)所示的关系：
n0＝nz*ny*x0+nz*y0+z0 (10)
fR1如式(11)所示：
Figure BDA0003917416640000093
式中，θr为参数，一般情况下取值0.2；
常数项θt0如公式(12)所示：
Figure BDA0003917416640000094
式中，常数p为目标电容率初值的数量；
S35、约束项R2为一规模为1*(nx*ny*nz)的列向量，约束项R2如公式(13)所示：
Figure BDA0003917416640000095
式中，λR2为约束项R2的系数，
Figure BDA0003917416640000104
为拉普拉斯算子，u为一规模为1*(nx*ny*nz)的列向量，u的每一项通式如公式(14)所示：
Figure BDA0003917416640000101
在上述步骤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的规模为nx*ny*nz，nx、ny、nz表示x，y，z方向上像素的个数，那么整个空间的像素个数即为nx*ny*nz。
在上述步骤S34中，目标电导率初值的数量与目标电容率初值的数量相同，因此常数p也可以为目标电导率初值的数量，例如，在步骤S31和步骤S32中筛选出1000个合理的目标电容率初值和目标电导率初值，即目标电容率初值和目标电导率初值的数量均为1000，因此常数p为1000。
在上述步骤S4中，基于发射场和发射场的点数据信息，将原始麦克斯韦方程和亥姆霍兹方程推导得到的高阶非线性非齐次偏微分方程通过改良的牛顿迭代法实现瞬时线化得到迭代公式，并根据物理模型和现实中电容率初值ε0和电导率初值σ0的情况加入约束项K、R1与R2，得到迭代方程
Figure BDA0003917416640000102
进一步地，在迭代方程中，f(u)如公式(2)所示：
f(u)＝Au+Cu2 (2)
式中，u2为一列向量，u2的每一项的通式为
Figure BDA0003917416640000103
A为由发射场和发射场的点数据信息确定的矩阵，如公式(3)所示：
Figure BDA0003917416640000111
式中，变量
Figure BDA0003917416640000116
如公式(4)所示，变量γ如公式(5)所示：
Figure BDA0003917416640000112
Figure BDA0003917416640000113
C为由发射场和发射场的点数据信息确定的矩阵，如公式(6)所示：
Figure BDA0003917416640000114
进一步地，在迭代方程中，阻尼因子向量α的计算方式如下：
构建一个nx*ny*nz的矩阵Q，Q每一点的值如公式(7)所示：
Figure BDA0003917416640000115
式中，σm和εm均为参数，(xi,yi,zi)与角标i一一对应，i表示第i个点；
设一矩阵V，对于V中每一点v,均有在Q中对应的以v坐标为中心的m*m矩阵T，v的值如公式(23)所示：
v＝min(T) (8)
将矩阵V一一映射到阻尼因子向量α中即得到阻尼因子向量α。
矩阵V的元素v对应空间中的每一个点，矩阵T相当于一个扫描窗口，其行数和列数m可以根据实际情况确定，优选地，m为5。而对于5*5矩阵T，每一个元素v，在矩阵Q中以对应的v坐标的点为中心点，以中心点和中心点各方向上的两个点构建一个5*5矩阵T，在构建的5*5矩阵T中找到最小的值作为元素v的值，以此不断得到每一个元素v的值，进而得到矩阵V。并且由上述过程可知，矩阵V在每次迭代时会改变，因此每次迭代后需要重新计算矩阵V，以重新计算阻尼因子向量α。
在上述步骤S5中，在迭代法中有一个合理的初值以加快收敛速度，因此通过根据空间中每一处的电容率初值和电导率初值，确定迭代方程的初值可以提高收敛速度。而本发明中磁共振使用的设备是正交线圈，故发射场B在z方向上分量可忽略不计，故所涉及原始方程如公式(23)所示：
Figure BDA0003917416640000121
式中，i表示虚数，ω为核磁共振设备对应的拉莫尔进动频率，μ0为磁导率，B′x为yz平面的发射场B沿x轴的偏导数，B′y为xz平面的发射场B沿y轴的偏导数，B′z为xy平面的发射场B沿z轴的偏导数；
根据Helmholtz EPT方法，假设电特性参数在空间中分布均匀，即空间中每一点处的电特性参数沿各方向导数均为0，使得公式(23)可以简化为如公式(24)所示：
Figure BDA0003917416640000122
由公式(24)推导可得到公式(25)：
Figure BDA0003917416640000123
进而得到迭代方程所需要的初值u0＝σ0+iωε0，构造的初值u0时，是从一个由空间中每一处的电容率初值和电导率初值构成的矩阵变换为一个向量u0，具有一一对应的关系，将初值u0代入迭代方程中进行迭代计算，经过若干次迭代后得到空间中每一点处的电容率和电导率。
进一步地，基于空间中每一处的电容率初值和电导率初值，确定迭代方程的初值，将初值代入迭代方程中进行迭代计算，经过若干次迭代后得到空间中每一点处的电容率和电导率的步骤包括：
S51、对迭代电容率初值和迭代电导率初值进行预处理，基于预处理后的迭代电容率初值和迭代电导率初值，根据公式(15)计算迭代方程的初值u0：
u0＝σ0+iωε0 (15)
式中，σ0为迭代电导率初，ε0为迭代电容率初值；
S52、将初值u0代入迭代方程中，得到新的初值，判断新的初值中的各项平均值是否小于预设的迭代阈值；
S53、若否，则重复步骤S52；
S54、若是，则根据迭代后新的初值、目标电容率初值和目标电导率初值得到空间中每一点处的电容率和电导率。
在上述步骤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，根据事先记录的坐标信息还原成矩阵形式，即可得到空间中每一点处的电容率值ε和电导率值σ。
在上述步骤S6中，空间中每一点处的电容率值ε和电导率值σ即为图像上每一点的信息，即可输出电特性断层成像图。
以下是本发明基于瞬时线化的电特性断层成像方法(以下简称本发明)与Helmholtz EPT和改良的crEPT方法的结果对比：
参考图为Duke人脑肿瘤仿真模型，评估指标选择峰值信噪比以及结构相似性来比较在本例中本发明与Helmholtz EPT和改良的crEPT方法的结果。峰值信噪比与结构相似性的值越高，说明重建出的电特性断层成像图与参考图像越接近，表1是本发明与HelmholtzEPT和改良的crEPT方法得出的结果的量化指标。
表1
方法 名称 峰值信噪比 结构相似性
Helmholtz EPT u 8.6151 0.6825
改良crEPT u 15.5295 0.8178
本发明 u 18.2423 0.8684
从表1可见，在本例中，无论是峰值信噪比还是结构相似性本发明均优于Helmholtz EPT和改良的crEPT方法。
另外结合请结合参阅图2和图3，图2为本发明基于瞬时线化的电特性断层成像方法与Helmholtz EPT方法、改良crEPT方法在Duke人脑肿瘤仿真模型上的重建电导率(上)与电容率(下)图像的结果对比，图3为本发明基于瞬时线化的电特性断层成像方法与Helmholtz EPT方法，改良crEPT方法在另一例Duke人脑肿瘤仿真模型上的重建电导率(上)与电容率(下)图像的结果。由图2和图3可知，本发明图视觉上更接近于参考图。因此，本发明得到的电特性断层成像图的效果要显著好于Helmholtz EPT和改良的crEPT方法。
请结合参阅图4，图4为本发明基于瞬时线化的电特性断层成像系统的模块示意图。本发明还提供了一种基于瞬时线化的电特性断层成像系统，包括：
测量模块1，用于测量磁共振射频场，基于射频场计算得到发射场；
第一计算模块2，用于利用数值方法计算得到发射场的点数据信息，点数据信息包括发射场在空间中每一点处沿x轴、y轴以及z轴方向的一阶导数和二阶导数；
第二计算模块3，用于基于发射场和发射场的点数据信息，根据麦克斯韦方程组，计算得到空间中每一点处的电容率初值和电导率初值，并由空间中每一点处的电容率初值和电导率初值确定约束项K、R1与R2；
得到模块4，用于基于发射场、发射场的点数据信息以及约束项K、R1与R2，根据麦克斯韦方程组结合牛顿迭代法，得到迭代方程，迭代方程如公式(1)所示：
Figure BDA0003917416640000142
式中，un表示第n代的u，即第n次迭代后得到的迭代值，f(u)为由发射场B的点数据信息确定的函数，α为阻尼因子向量；
迭代模块5，用于基于空间中每一处的电容率初值和电导率初值，确定迭代方程的初值，将初值代入迭代方程中进行迭代计算，经过若干次迭代后得到空间中每一点处的电容率和电导率；
输出模块6，用于基于空间中每一点处的电容率和电导率输出电特性断层成像图。
在上述测量模块1中，用核磁共振设备配合正交线圈以及特定B1扫描序列可以得到正交线圈测量空间的射频场B，即得到磁共振射频场B。然后基于得到的射频场B，计算出发射场B。
进一步地，测量模块1包括：
测量子模块，测量磁共振射频场沿空间坐标系x轴、y轴以及z轴方向的分量Bx、By、Bz；
第一计算子模块，用于根据公式(22)计算发射场B：
Figure BDA0003917416640000141
式中，i为复数单位。
在上述测量子模块至第一计算子模块中，由仪器测得的射频场B唯一矢量，其值为(Bx,By,Bz)，其中，Bx、By、Bz分别为射频场B沿空间坐标系x轴、y轴以及z轴方向的分量。而发射场B由公式
Figure BDA0003917416640000151
计算。
数值方法可以采用有限差分方法，基于发射场B，第一计算模块2通过有限差分方法计算得到发射场在空间中每一点处沿x轴、y轴以及z轴方向的一阶导数和二阶导数。
进一步地，第一计算模块2包括：
第二计算子模块，用于根据公式(16)计算发射场在空间中每一点处沿x轴方向的一阶导数
Figure BDA0003917416640000152
Figure BDA0003917416640000153
第三计算子模块，用于根据公式(17)计算发射场在空间中每一点处沿y轴方向的一阶导数
Figure BDA0003917416640000154
Figure BDA0003917416640000155
第四计算子模块，用于根据公式(18)计算发射场在空间中每一点处沿z轴方向的一阶导数
Figure BDA0003917416640000156
Figure BDA0003917416640000157
第五计算子模块，用于根据公式(19)计算发射场在空间中每一点处沿x轴方向的二阶导数
Figure BDA0003917416640000158
Figure BDA0003917416640000159
第六计算子模块，用于根据公式(20)计算发射场在空间中每一点处沿y轴方向的二阶导数
Figure BDA00039174166400001510
Figure BDA00039174166400001511
第七计算子模块，用于根据公式(21)计算发射场在空间中每一点处沿z轴方向的二阶导数
Figure BDA00039174166400001512
Figure BDA00039174166400001513
其中，dx、dy、dz分别为空间离散化后的单位网格在x轴、y轴以及z轴方向上的长度，i为角标，表示空间中特定一点。
第二计算模块3基于发射场和发射场的点数据信息，根据MREPT(magneticresonance electrical properties，人体组织电特性磁共振断层成像)的原理，计算得到空间中每一点处的电容率初值ε0和电导率初值σ0，然后根据预设的处理规则对计算得到空间中每一点处的电容率初值ε0和电导率初值σ0进行处理，确定约束项K、R1与R2。
进一步地，第二计算模块3包括：
判断子模块，用于对于空间中每一点处的电容率初值和电导率初值，判断电容率初值和电导率初值沿各方向导数是否为0；
第一标记子模块，用于若判断子模块的判断为是，则将电容率初值和电导率初值记为目标电容率初值和目标电导率初值，目标电容率初值和目标电导率初值不参与迭代，并根据目标电容率初值和目标电导率初值构造约束项K；
第二标记子模块，用于若判断子模块的判断为否，则将电容率初值和电导率初值记为迭代电容率初值和迭代电导率初值；
第一约束项模块，用于约束项R1为规模为1*(nx*ny*nz)的列向量，约束项R1中每一项的通式如公式(9)所示：
Figure BDA0003917416640000161
式中，λR1为约束项R1的系数，
Figure BDA0003917416640000162
ω为核磁共振设备对应的拉莫尔进动频率，n0与(x0,y0,z0)满足如公式(10)所示的关系：
n0＝nz*ny*x0+nz*y0+z0 (10)
fR1如式(11)所示：
Figure BDA0003917416640000163
式中，θr为参数，一般情况下取值0.2；
常数项θt0如公式(12)所示：
Figure BDA0003917416640000164
式中，常数p为目标电容率初值的数量；
第一约束项模块，用于约束项R2为一规模为1*(nx*ny*nz)的列向量，约束项R2如公式(13)所示：
Figure BDA0003917416640000165
式中，λR2为约束项R2的系数，
Figure BDA0003917416640000167
为拉普拉斯算子，u为一规模为1*(nx*ny*nz)的列向量，u的每一项通式如公式(14)所示：
Figure BDA0003917416640000166
参考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的规模为nx*ny*nz，nx、ny、nz表示x，y，z方向上像素的个数，那么整个空间的像素个数即为nx*ny*nz。
在上述第一约束项模块中，目标电导率初值的数量与目标电容率初值的数量相同，因此常数p也可以为目标电导率初值的数量，例如，由判断子模块和第一标记子模块筛选出1000个合理的目标电容率初值和目标电导率初值，即目标电容率初值和目标电导率初值的数量均为1000，因此常数p为1000。
得到模块4基于发射场和发射场的点数据信息，将原始麦克斯韦方程和亥姆霍兹方程推导得到的高阶非线性非齐次偏微分方程通过改良的牛顿迭代法实现瞬时线化得到迭代公式，并根据物理模型和现实中电容率初值ε0和电导率初值σ0的情况加入约束项K、R1与R2，得到迭代方程
Figure BDA0003917416640000171
进一步地，在迭代方程中，f(u)如公式(2)所示：
f(u)＝Au+Cu2 (2)
式中，u2为一列向量，u2的每一项的通式为
Figure BDA0003917416640000172
A为由发射场和发射场的点数据信息确定的矩阵，如公式(3)所示：
Figure BDA0003917416640000181
式中，变量
Figure BDA0003917416640000186
如公式(4)所示，变量γ如公式(5)所示：
Figure BDA0003917416640000182
Figure BDA0003917416640000183
C为由发射场和发射场的点数据信息确定的矩阵，如公式(6)所示：
Figure BDA0003917416640000184
进一步地，在迭代方程中，阻尼因子向量α的计算方式如下：
构建一个nx*ny*nz的矩阵Q，Q每一点的值如公式(7)所示：
Figure BDA0003917416640000185
式中，σm和εm均为参数，(xi,yi,zi)与角标i一一对应，i表示第i个点；
设一矩阵V，对于V中每一点v,均有在Q中对应的以v坐标为中心的m*m矩阵T，v的值如公式(23)所示：
v＝min(T) (8)
将矩阵V一一映射到阻尼因子向量α中即得到阻尼因子向量α。
矩阵V的元素v对应空间中的每一个点，矩阵T相当于一个扫描窗口，其行数和列数m可以根据实际情况确定，优选地，m为5。而对于5*5矩阵T，每一个元素v，在矩阵Q中以对应的v坐标的点为中心点，以中心点和中心点各方向上的两个点构建一个5*5矩阵T，在构建的5*5矩阵T中找到最小的值作为元素v的值，以此不断得到每一个元素v的值，进而得到矩阵V。并且由上述过程可知，矩阵V在每次迭代时会改变，因此每次迭代后需要重新计算矩阵V，以重新计算阻尼因子向量α。
在上述迭代模块5中，在迭代法中有一个合理的初值以加快收敛速度，因此通过根据空间中每一处的电容率初值和电导率初值，确定迭代方程的初值可以提高收敛速度。而本发明中磁共振使用的设备是正交线圈，故发射场B在z方向上分量可忽略不计，故所涉及原始方程如公式(23)所示：
Figure BDA0003917416640000191
式中，i表示虚数，ω为核磁共振设备对应的拉莫尔进动频率，μ0为磁导率，B′x为yz平面的发射场B沿x轴的偏导数，B′y为xz平面的发射场B沿y轴的偏导数，B′z为xy平面的发射场B沿z轴的偏导数；
根据Helmholtz EPT方法，假设电特性参数在空间中分布均匀，即空间中每一点处的电特性参数沿各方向导数均为0，使得公式(23)可以简化为如公式(24)所示：
Figure BDA0003917416640000192
由公式(24)推导可得到公式(25)：
Figure BDA0003917416640000193
进而得到迭代方程所需要的初值u0＝σ0+iωε0，构造的初值u0时，是从一个由空间中每一处的电容率初值和电导率初值构成的矩阵变换为一个向量u0，具有一一对应的关系，将初值u0代入迭代方程中进行迭代计算，经过若干次迭代后得到空间中每一点处的电容率和电导率。
进一步地，迭代模块5包括：
处理子模块，用于对迭代电容率初值和迭代电导率初值进行预处理，基于预处理后的迭代电容率初值和迭代电导率初值，根据公式(15)计算迭代方程的初值u0：
u0＝σ0+iωε0 (15)
式中，σ0为迭代电导率初，ε0为迭代电容率初值；
代入模块，用于将初值u0代入迭代方程中，得到新的初值，判断新的初值中的各项平均值是否小于预设的迭代阈值；
重复子模块，用于若否，则重复步骤S52；
得到子模块，用于若是，则根据迭代后新的初值、目标电容率初值和目标电导率初值得到空间中每一点处的电容率和电导率。
在上述处理子模块至得到子模块中，由判断子模块、第一标记子模块和第二标记子模块筛选出需要进行迭代的迭代电容率初值εd0和电导率初值σd0，根据公式(15)构成初值u0，将初值u0代入迭代方程中，得到新的初值u1，新的初值u1中的各项平均值是否小于预设的迭代阈值，判断f(u1)的均值是否小迭代阈值，若否，则根据新的初值u1代入迭代方程，得到新的初值u2，以此不断进行迭代，直到某一代的初值ui中的各项平均值小于预设的迭代阈值，则结束迭代，得到ui。由于构件初值u0时，是从一个矩阵变换为一个向量，具有一一对应的关系，因此可以事先记录每个迭代电容率初值εd0和电导率初值σd0的映射信息，及空间中对应的坐标，同样事先记录由步骤S31至步骤S33筛选出来不需要参与迭代的目标电容率初值εt0和目标电导率初值σt0的映射信息，在计算得到ui后，基于ui、目标电容率初值εt0和目标电导率初值σt0，根据事先记录的坐标信息还原成矩阵形式，即可得到空间中每一点处的电容率值ε和电导率值σ。
在上述输出模块6中，空间中每一点处的电容率值ε和电导率值σ即为图像上每一点的信息，即可输出电特性断层成像图。
请结合参阅图5，图5为本发明电子设备一实施例的结构示意框图。本发明一实施例还提出一种电子设备1001，包括存储器1003和处理器1002，存储器1003存储有计算机程序1004，处理器1002执行计算机程序1004时实现上述任一项基于瞬时线化的电特性断层成像方法的步骤，包括：S1、测量磁共振射频场，基于射频场计算得到发射场；S2、利用数值方法计算得到发射场的点数据信息，点数据信息包括发射场在空间中每一点处沿x轴、y轴以及z轴方向的一阶导数和二阶导数；S3、基于发射场和发射场的点数据信息，根据麦克斯韦方程组，计算得到空间中每一点处的电容率初值和电导率初值，并由空间中每一点处的电容率初值和电导率初值确定约束项K、R1与R2；S4、基于发射场、发射场的点数据信息以及约束项K、R1与R2，根据麦克斯韦方程组结合牛顿迭代法，得到迭代方程
Figure BDA0003917416640000201
Figure BDA0003917416640000202
S5、基于空间中每一处的电容率初值和电导率初值，确定迭代方程的初值，将初值代入迭代方程中进行迭代计算，经过若干次迭代后得到空间中每一点处的电容率和电导率；S6、基于空间中每一点处的电容率和电导率输出电特性断层成像图。
请结合参阅图6，图6为本发明计算机可读存储介质一实施例的结构示意框图。本发明一实施例还提供一种计算机可读存储介质2001，其上存储有计算机程序1004，计算机程序1004被处理器1002执行时实现上述任一项基于瞬时线化的电特性断层成像方法的步骤，包括：S1、测量磁共振射频场，基于射频场计算得到发射场；S2、利用数值方法计算得到发射场的点数据信息，点数据信息包括发射场在空间中每一点处沿x轴、y轴以及z轴方向的一阶导数和二阶导数；S3、基于发射场和发射场的点数据信息，根据麦克斯韦方程组，计算得到空间中每一点处的电容率初值和电导率初值，并由空间中每一点处的电容率初值和电导率初值确定约束项K、R1与R2；S4、基于发射场、发射场的点数据信息以及约束项K、R1与R2，根据麦克斯韦方程组结合牛顿迭代法，得到迭代方程
Figure BDA0003917416640000211
Figure BDA0003917416640000212
S5、基于空间中每一处的电容率初值和电导率初值，确定迭代方程的初值，将初值代入迭代方程中进行迭代计算，经过若干次迭代后得到空间中每一点处的电容率和电导率；S6、基于空间中每一点处的电容率和电导率输出电特性断层成像图。
相比于现有技术，本发明的有益效果为：将原始麦克斯韦方程和亥姆霍兹方程推导得到的高阶非线性非齐次偏微分方程通过改良的牛顿迭代法实现瞬时线化得到迭代方程，并根据物理模型和现实情况加入约束项，构成完备的具有初值的迭代方程，得到一个精确度更高，更符合真实状况的电特性分布解，从而能够在尽量少简化或者更改问题模型的前提下，更加准确地得到电特性断层成像。
本领域普通技术人员可以理解实现上述实施例方法中的全部或部分流程，是可以通过计算机程序来指令相关的硬件来完成，所述的计算机程序可存储于一非易失性计算机可读取存储介质中，该计算机程序在执行时，可包括如上述各方法的实施例的流程。其中，本申请所提供的和实施例中所使用的对存储器、存储、数据库或其它介质的任何引用，均可包括非易失性和/或易失性存储器。非易失性存储器可以包括只读存储器(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)等。
需要说明的是，在本文中，术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性的包含，从而使得包括一系列要素的过程、装置、物品或者方法不仅包括那些要素，而且还包括没有明确列出的其他要素，或者是还包括为这种过程、装置、物品或者方法所固有的要素。在没有更多限制的情况下，由语句“包括一个……”限定的要素，并不排除在包括该要素的过程、装置、物品或者方法中还存在另外的相同要素。
本发明并不局限于上述实施方式，如果对本发明的各种改动或变形不脱离本发明的精神和范围，倘若这些改动和变形属于本发明的权利要求和等同技术范围之内，则本发明也意图包含这些改动和变形。

Claims (10)

1.一种基于瞬时线化的电特性断层成像方法，其特征在于，包括以下步骤：
S1、测量磁共振射频场，基于射频场计算得到发射场；
S2、利用数值方法计算得到发射场的点数据信息，所述点数据信息包括发射场在空间中每一点处沿x轴、y轴以及z轴方向的一阶导数和二阶导数；
S3、基于发射场和发射场的点数据信息，根据麦克斯韦方程组，计算得到空间中每一点处的电容率初值和电导率初值，并由空间中每一点处的电容率初值和电导率初值确定约束项K、R1与R2；
S4、基于发射场、发射场的点数据信息以及约束项K、R1与R2，根据麦克斯韦方程组结合牛顿迭代法，得到迭代方程，所述迭代方程如公式(1)所示：
Figure FDA0003917416630000011
式中，un表示第n代的u，f(u)为由发射场B的点数据信息确定的函数，α为阻尼因子向量；
S5、基于空间中每一处的电容率初值和电导率初值，确定迭代方程的初值，将初值代入迭代方程中进行迭代计算，经过若干次迭代后得到空间中每一点处的电容率和电导率；
S6、基于空间中每一点处的电容率和电导率输出电特性断层成像图。
2.根据权利要求1所述的基于瞬时线化的电特性断层成像方法，其特征在于，在迭代方程中，f(u)如公式(2)所示：
f(u)＝Au+Cu2 (2)
式中，u2为一列向量，u2的每一项的通式为
Figure FDA0003917416630000012
A为矩阵，如公式(3)所示：
Figure FDA0003917416630000013
式中，变量
Figure FDA0003917416630000021
如公式(4)所示，变量γ如公式(5)所示：
Figure FDA0003917416630000022
Figure FDA0003917416630000023
C为矩阵，如公式(6)所示：
Figure FDA0003917416630000024
3.根据权利要求1所述的基于瞬时线化的电特性断层成像方法，其特征在于，在迭代方程中，阻尼因子向量α的计算方式如下：
构建一个nx*ny*nz的矩阵Q，Q每一点的值如公式(7)所示：
Figure FDA0003917416630000025
式中，σm和εm均为参数；
设一矩阵V，对于V中每一点v,均有在Q中对应的以v坐标为中心的m*m矩阵T，v的值如公式(23)所示：
v＝min(T) (8)
将矩阵V一一映射到阻尼因子向量α中即得到阻尼因子向量α。
4.根据权利要求1所述的基于瞬时线化的电特性断层成像方法，其特征在于，所述由空间中每一点处的电容率初值和电导率初值确定约束项K、R1与R2的步骤包括：
对于空间中每一点处的电容率初值和电导率初值，判断电容率初值和电导率初值沿各方向导数是否为0；
若是，则将电容率初值和电导率初值记为目标电容率初值和目标电导率初值，目标电容率初值和目标电导率初值不参与迭代，并根据目标电容率初值和目标电导率初值构造约束项K；
若否，则将电容率初值和电导率初值记为迭代电容率初值和迭代电导率初值；
约束项R1为规模为1*(nx*ny*nz)的列向量，约束项R1中每一项的通式如公式(9)所示：
Figure FDA0003917416630000031
式中，λR1为约束项R1的系数，
Figure FDA0003917416630000032
ω为核磁共振设备对应的拉莫尔进动频率，n0与(x0,y0,z0)满足如公式(10)所示的关系：
n0＝nz*ny*x0+nz*y0+z0 (10)
fR1如式(11)所示：
Figure FDA0003917416630000033
式中，θr为参数；
常数项θt0如公式(12)所示：
Figure FDA0003917416630000034
式中，常数p为目标电容率初值的数量；
约束项R2为一规模为1*(nx*ny*nz)的列向量，约束项R2如公式(13)所示：
Figure FDA0003917416630000035
式中，λR2为约束项R2的系数，
Figure FDA0003917416630000036
为拉普拉斯算子，u为一规模为1*(nx*ny*nz)的列向量，u的每一项通式如公式(14)所示：
Figure FDA0003917416630000037
5.根据权利要求4所述的基于瞬时线化的电特性断层成像方法，其特征在于，所述基于空间中每一处的电容率初值和电导率初值，确定迭代方程的初值，将初值代入迭代方程中进行迭代计算，经过若干次迭代后得到空间中每一点处的电容率和电导率的步骤包括：
S51、对迭代电容率初值和迭代电导率初值进行预处理，基于预处理后的迭代电容率初值和迭代电导率初值，根据公式(15)计算迭代方程的初值u0：
u0＝σ0+iωε0 (15)
式中，σ0为迭代电导率初，ε0为迭代电容率初值；
S52、将初值u0代入迭代方程中，得到新的初值，判断新的初值中的各项平均值是否小于预设的迭代阈值；
S53、若否，则重复步骤S52；
S54、若是，则根据迭代后新的初值、目标电容率初值和目标电导率初值得到空间中每一点处的电容率和电导率。
6.根据权利要求1所述的基于瞬时线化的电特性断层成像方法，其特征在于，所述利用数值方法计算得到发射场的点数据信息的步骤包括：
根据公式(16)计算发射场在空间中每一点处沿x轴方向的一阶导数
Figure FDA0003917416630000041
Figure FDA0003917416630000042
根据公式(17)计算发射场在空间中每一点处沿y轴方向的一阶导数
Figure FDA0003917416630000043
Figure FDA0003917416630000044
根据公式(18)计算发射场在空间中每一点处沿z轴方向的一阶导数
Figure FDA0003917416630000045
Figure FDA0003917416630000046
根据公式(19)计算发射场在空间中每一点处沿x轴方向的二阶导数
Figure FDA0003917416630000047
Figure FDA0003917416630000048
根据公式(20)计算发射场在空间中每一点处沿y轴方向的二阶导数
Figure FDA0003917416630000049
Figure FDA00039174166300000410
根据公式(21)计算发射场在空间中每一点处沿z轴方向的二阶导数
Figure FDA00039174166300000411
Figure FDA00039174166300000412
其中，dx、dy、dz分别为空间离散化后的单位网格在x轴、y轴以及z轴方向上的长度。
7.根据权利要求1所述的基于瞬时线化的电特性断层成像方法，其特征在于，所述测量磁共振射频场，基于射频场计算得到发射场的步骤包括：
测量磁共振射频场沿空间坐标系x轴、y轴以及z轴方向的分量Bx、By、Bz；
根据公式(22)计算发射场B：
Figure FDA0003917416630000051
式中，i为复数单位。
8.一种基于瞬时线化的电特性断层成像系统，其特征在于，包括：
测量模块，用于测量磁共振射频场，基于射频场计算得到发射场；
第一计算模块，用于利用数值方法计算得到发射场的点数据信息，所述点数据信息包括发射场在空间中每一点处沿x轴、y轴以及z轴方向的一阶导数和二阶导数；
第二计算模块，用于基于发射场和发射场的点数据信息，根据麦克斯韦方程组，计算得到空间中每一点处的电容率初值和电导率初值，并由空间中每一点处的电容率初值和电导率初值确定约束项K、R1与R2；
得到模块，用于基于发射场、发射场的点数据信息以及约束项K、R1与R2，根据麦克斯韦方程组结合牛顿迭代法，得到迭代方程，所述迭代方程如公式(1)所示：
Figure FDA0003917416630000052
式中，un表示第n代的u，f(u)为由发射场B的点数据信息确定的函数，α为阻尼因子向量；
迭代模块，用于基于空间中每一处的电容率初值和电导率初值，确定迭代方程的初值，将初值代入迭代方程中进行迭代计算，经过若干次迭代后得到空间中每一点处的电容率和电导率；
输出模块，用于基于空间中每一点处的电容率和电导率输出电特性断层成像图。
9.一种电子设备，包括存储器和处理器，所述存储器存储有计算机程序，其特征在于，所述处理器执行所述计算机程序时实现权力要求1-7任一项所述方法的步骤。
10.一种计算机可读存储介质，其上存储有计算机程序，其特征在于，所述计算机程序被处理器执行时实现权利要求1-7任一项所述的方法的步骤。
CN202211343568.4A 2022-10-31 2022-10-31 基于瞬时线化的电特性断层成像方法、系统、设备及介质 Pending CN115808650A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202211343568.4A CN115808650A (zh) 2022-10-31 2022-10-31 基于瞬时线化的电特性断层成像方法、系统、设备及介质

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202211343568.4A CN115808650A (zh) 2022-10-31 2022-10-31 基于瞬时线化的电特性断层成像方法、系统、设备及介质

Publications (1)

Publication Number Publication Date
CN115808650A true CN115808650A (zh) 2023-03-17

Family

ID=85482867

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202211343568.4A Pending CN115808650A (zh) 2022-10-31 2022-10-31 基于瞬时线化的电特性断层成像方法、系统、设备及介质

Country Status (1)

Country Link
CN (1) CN115808650A (zh)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116544931A (zh) * 2023-06-27 2023-08-04 北京理工大学 基于集成片段变换和时间卷积网络的电力负荷分布预测方法
CN117113794A (zh) * 2023-10-23 2023-11-24 之江实验室 磁约束带电粒子成像系统中反角度准直器的设计方法

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116544931A (zh) * 2023-06-27 2023-08-04 北京理工大学 基于集成片段变换和时间卷积网络的电力负荷分布预测方法
CN116544931B (zh) * 2023-06-27 2023-12-01 北京理工大学 基于集成片段变换和时间卷积网络的电力负荷分布预测方法
CN117113794A (zh) * 2023-10-23 2023-11-24 之江实验室 磁约束带电粒子成像系统中反角度准直器的设计方法
CN117113794B (zh) * 2023-10-23 2024-01-26 之江实验室 磁约束带电粒子成像系统中反角度准直器的设计方法

Similar Documents

Publication Publication Date Title
CN115808650A (zh) 基于瞬时线化的电特性断层成像方法、系统、设备及介质
US8103337B2 (en) Weighted gradient method and system for diagnosing disease
Gao et al. Sensitivity of the distorted born iterative method to the initial guess in microwave breast imaging
Joachimowicz et al. Convergence and stability assessment of Newton-Kantorovich reconstruction algorithms for microwave tomography
CN101331516B (zh) 用于多次迭代算法的高级收敛
EP1520241B1 (en) Method and system for displaying confidence intervals for source reconstruction
CN114270397A (zh) 使用电特性断层成像确定流体和组织体积估计的系统和方法
CN104851080A (zh) 一种基于tv的三维pet图像重建方法
Huttinga et al. Gaussian Processes for real-time 3D motion and uncertainty estimation during MR-guided radiotherapy
Singh et al. 3-D EM modeling of medical microwave imaging scenarios with controllable accuracy
AU2004224834A1 (en) Weighted gradient method and system for diagnosing disease
EP3833246A1 (en) System, method, and computer-accessible medium for non-invasive temperature estimation
Lacik et al. Rat head phantom for testing of electroencephalogram source localization techniques
CN110720914A (zh) 基于稀疏采样的全息磁感应胸腔成像方法及成像系统
Plis et al. Probabilistic forward model for electroencephalography source analysis
JP3757276B2 (ja) 磁気共鳴イメージング装置の画像情報処理方法、及び磁気共鳴イメージング装置
Meeson et al. The dependence of EIT images on the assumed initial conductivity distribution: a study of pelvic imaging
Hashemzadeh et al. A hybrid analytical–numerical algorithm for determining the neuronal current via electroencephalography
DE112021006710T5 (de) Qualitätsmass für eine mappingfunktion
KR20210068189A (ko) 의료 영상을 이용한 병변 여부 판단 방법
Caeiros et al. A new image reconstruction algorithm for real-time monitoring of conductivity and permeability changes in magnetic induction tomography
US20040243019A1 (en) Weighted gradient method and system for diagnosing disease
CN109859173B (zh) 一种基于电磁逆散射的早期乳腺癌检测医学成像方法
Golyshev et al. Analysis of methods for localization magnetic field sources in biomagnetic studies
Zaravi et al. Expanding 2D block method in two direction by a new formula in EIT

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination