一种基于物理神经网络的电特性断层成像方法

1.本发明涉及核磁共振成像技术领域，特别涉及一种基于物理神经网络的电特性断层成像方法。


背景技术：

2.高场核磁共振成像系统随着主场场强的提高，射频场的频率变高，其波长会变短到跟被扫描物体在同一个水平。这会引起射频场和被扫描物体的相互作用增强。从而可以利用这种相互作用实现一些新的成像功能。基于高场磁共振的介电特性断层成像(mr-ept)就是一种新的核磁共振量化成像方法。
3.生物组织介电特性是生物组织作为一种物质对电磁场的响应特性。介电特性不仅可以作为组织对比度显示人体内不同组织的结构图像，还可以作为一种定量指标，用于早期癌症的检测。除此之外，介电特性会影响脑内电信号的传播，并与脑组织的类型、活动、功能、衰老、病理改变等信息高度，相关因此mr-ept有望为研究脑认知与脑疾病提供全新的研究工具。人体内介电特性分布还是计算特定吸收率(sar)的基础，用于评估任何暴露于射频/微波的射频安全。超高场磁共振成像系统面临的一个很大的安全问题，就是生物组织在强磁场中的热效应。对于人体组织的介电特性的准确成像能够有助于sar的研究，推动超高场磁共振系统的进一步发展。
4.mr-ept一般分为两个步骤，第一步是射频场幅值和相位的采集。第二步是以麦克斯韦方程为基础，得到mr-ept的核心方程，然后构造逆问题的数学模型，并利用已采集的射频场数据作为逆问题的输入，来重建电特性的分布。其中第二步是电特性成像算法的核心，尤其是如何对mr-ept的核心方程进行重构来形成数学模型，并根据所构造的数学模型对电特性分布进行精确的重建。
5.目前，用于电特性重建的主流算法，都是利用一些假设来对mr-ept的核心方程进行简化，从而形成简单的求解方程。有一种假设是认为被测物体内部的电特性分布是均匀的，这种方法被称为helmholtz ept。虽然这种假设可以极大的简化求解模型，可以快速的估算出电特性分布，但是利用这种均匀性假设有着明显的缺点，因为被测物体内部的电特性一般都会具有较复杂的分布，而不是均匀一块。这会使得均匀假设失效，从而导致严重的计算误差，特别是在组织边缘部分产生严重的伪影。同时，亥姆霍兹方程本身是一个高阶非线性非齐次偏微分方程，并且射频场和电特性均以复数形式体现在方程当中。所以在不利用假设来简化的情况下，很难直接对方程求解。
6.因此，针对现有技术不足，提供一种纳米纤维表面多孔骨修复支架及其制备方法以解决现有技术不足甚为必要。


技术实现要素：

7.本发明的目的在于避免现有技术的不足之处而提供一种基于物理神经网络的电特性断层成像方法。该基于物理神经网络的电特性断层成像方法使用神经网络构建了从发
射场以及其空间导数到电特性之间的映射关系，得到结果的精度更高。
8.本发明的上述目的通过以下技术措施实现：
9.提供一种基于物理神经网络的电特性断层成像方法，包括如下步骤：
10.步骤(1)、根据射频场b，得到发射场b；
11.步骤(2)、利用数值方法求解所述步骤(1)得到的发射场b中每点在空间x轴、y轴和z轴上的一阶导数及二阶导数；
12.步骤(3)、构建全连接神经网络，并对麦克斯韦方程进行组合变换得到损失函数；
13.步骤(4)、使用发射场b、步骤(2)得到的一阶导数和二阶导数构成的数据集及步骤(3)得到的损失函数训练步骤(3)得到的全连接神经网络，得到训练后全连接神经网络，所述训练后全连接神经网络以空间任意一点的发射场b以及通过步骤(2)求出在该点的一阶导数为输入，所述训练后全连接神经网络的输出为该点的电容率及电导率，从而得到空间所有点的电容率及电导率。
14.在所述步骤(1)中，射频场b为射频线圈产生的，且b＝(b
x
,by,bz)，其中b
x
为射频场在空间上x轴的分量，by为射频场在空间上y轴的分量，bz为射频场在空间上z轴的分量，根据射频场b得到发射场b，存在
15.所述步骤(3)中，损失函数的获得方法步骤具体为：
16.步骤(3.1)、对麦克斯韦方程组进行组合变换，得到式(1)
[0017][0018]
其中εc＝ωε-iσ为复介电常数，i为虚数单位，ε为电容率，σ为电导率，ω为核磁共振设备对应的拉莫尔进动频率，μ为磁导率，为nabla算子，为拉普拉斯算子；
[0019]
步骤(3.2)、令bz为0并对步骤(1.1)的式(1)进行分离和重组，得到mr-ept的核心方程式(2)；
[0020][0021]
其中b
x
为发射场b中在空间上x轴的一阶导数,by为发射场b中在空间上y轴的一阶导数,bz为发射场b中在空间上z轴的一阶导数,ε
c,x
为复介电常数在空间上x轴上的一阶导数，ε
c,y
为复介电常数在空间上y轴上的一阶导数，pde为偏微分方程符号；
[0022]
步骤(3.3)、根据式(2)得到损失函数。
[0023]
优选的，上述损失函数由式(3)表示，
[0024][0025]
其中由式(4)表示，
[0026][0027]
其中由式(5)表示，
[0028][0029]
其中，α为平衡和的权重且α》0，corr函数为计算电容率梯度和电导率梯度之间的相关性之间的互相关函数，由式(6)表示，
[0030][0031]
其中为x的均值，为y的均值，n为计算点的个数，k为第k个点。
[0032]
优选的，上述步骤(2)具体为利用数值方法求解发射场b中每点在空间三个方向的一阶导数b
x
、一阶导数by、一阶导数bz、二阶导数b
xx
、二阶导数b
xy
、二阶导数b
xz
、二阶导数b
yz
、二阶导数b
yy
和二阶导数b
zz
；
[0033]
一阶导数b
x
由式(7)表示，
[0034][0035]
一阶导数by由式(8)表示，
[0036][0037]
一阶导数bz由式(9)表示,
[0038][0039]
二阶导数b
xx
由式(10)表示；
[0040][0041]
二阶导数b
xy
由式(11)表示，
[0042][0043]
二阶导数b
xz
由式(12)表示，
[0044][0045]
二阶导数b
yz
由式(13)表示，
[0046][0047]
二阶导数b
yy
由式(14)表示，
[0048][0049]
二阶导数b
zz
由式(15)表示；
[0050][0051]
其中(xk,yk,zk)为第k个点的坐标，dx,dy和dz分别为在空间三个方向上离散化后的网格。
[0052]
优选的，上述全连接神经网络分由网络a和网络b两部分组成。
[0053]
优选的，上述网络a构造了f
ε
函数，所述网络b构造了f
σ
函数，所述网络a和所述网络b输入均为空间上任意一点的发射场b以及该点的一阶导数b
x
、一阶导数by和一阶导数bz，输出分别为该点的电容率和电导率。
[0054]
优选的，上述网络a的输入层的激活函数、所述网络a的中间层的激活函数、所述网络b的输入层的激活函数及所述网络b的输入层的中间层的激活函数均为tanh函数，
[0055]
且tanh函数由式(16)表示，
[0056][0057]
优选的，上述网络a的输出层的激活函数、所述网络b的输出层的激活函数均为sigmoid函数，且sigmoid函数由式(17)表示，
[0058][0059]
其中tanh函数将输入映射到(-1,1)之间，sigmoid函数将输入映射到(0,1)之间。
[0060]
当所述网络a和所述网络b输入均为空间上任意一点的发射场b及该点的一阶导数b
x
、一阶导数by和一阶导数bz时，该点的电容率由式(18)表示，该点的电导率由式(19)表示；
[0061]
ε＝(81*oa 1)ecꢀꢀ
式(18)；
[0062]
σ＝3*obꢀꢀ
式(19)；
[0063]
其中ε为电导率，σ为电容率，oa为网络a的输出，ob为网络b的输出。
[0064]
优选的，上述全连接神经网络的映射由式(20)和式(21)表示，
[0065]
ε(xk,yk,zk)
[0066]
＝f
ε
(b(xk,yk,zk),b
x
(xk,yk,zk),by(xk,yk,zk),bz(xk,yk,zk)；w
ε
)
[0067]
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
式(20)；
[0068]
σ(xk,yk,zk)
[0069]
＝f
σ
(b(xk,yk,zk),b
x
(xk,yk,zk),by(xk,yk,zk),bz(xk,yk,zk)；w
σ
)
[0070]
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
式(21)；
[0071]
其中w
ε
和w
σ
均为神经网络的参数，ε(xk,yk,zk)为空间上坐标为(xk,yk,zk)上的电容率值，σ(xk,yk,zk)为空间上坐标为(xk,yk,zk)上的电导率值。
[0072]
优选的，上述步骤(4)具体为：
[0073]
步骤(4.1)、使用发射场b、步骤(2)得到的一阶导数和二阶导数构成的数据集全部输入步骤(3)得到的全连接神经网络，得到全连接神经网络的输出值，将全连接神经网络的输出值输入至步骤(3)得到的损失函数，令n＝1并进行前向和反向运算得到数据集的当前平均值平均值更新全连接神经网络，进入步骤(4.2)；
[0074]
步骤(4.2)、令n＝n 1，进入步骤(4.3)；
[0075]
步骤(4.3)、使用发射场b、步骤(2)得到的一阶导数和二阶导数构成的数据集全部输入更新后的全连接神经网络，得到全连接神经网络的输出值，进入步骤(4.4)；
[0076]
步骤(4.4)、将全连接神经网络的输出值输入至步骤(3)得到的损失函数，然后进行前向和反向运算得到数据集的平均值判断与迭代阈值的关系，当大于等迭代阈值则进入步骤(4.5)，否则进入步骤(4.6)；
[0077]
步骤(4.5)、将步骤(4.4)得到的平均值更新全连接神经网络，得到更新后的全连接神经网络并返回步骤(4.2)；
[0078]
步骤(4.6)、将当前的更新全连接神经网络为定义为训练后全连接神经网络。
[0079]
优选的，上述数值方法为有限差分方法；
[0080]
优选的，上述迭代阈值为0.001。
[0081]
优选的，上述网络a和所述网络b均由五层全连接层构成，且每层有256个神经元。
[0082]
优选的，上述全连接神经网络使用python3.6编写，并使用tensorflow深度学习库构建。
[0083]
一种基于物理神经网络的电特性断层成像方法，包括如下步骤：步骤(1)、根据射频场b，得到发射场b；步骤(2)、利用数值方法求解所述步骤(1)得到的发射场b中每点在空间三个方向的一阶导数及二阶导数；步骤(3)、构建全连接神经网络，并对麦克斯韦方程进行组合变换得到损失函数；步骤(4)、使用发射场b、步骤(2)得到的一阶导数和二阶导数构成的数据集、步骤(3)得到的损失函数训练步骤(3)得到的全连接神经网络，得到训练后全连接神经网络，所述训练后全连接神经网络以空间任意一点的发射场b以及通过步骤(2)求出在该点的一阶导数为输入，所述训练后全连接神经网络的输出为该点的电容率及电导率，从而得到空间所有点的电容率及电导率。本发明不需要对mr-ept的核心方程做简化来求解，而且在计算介电特性在空间的导数时不需要性离散化利用数值方法进行求解，避免了离散误差，因此本发明得到电容率及电导率的精度更高。该基于物理神经网络的电特性断层成像方法使用神经网络构建了从发射场以及其空间导数到电特性之间的映射关系，得到结果的精度更高。
附图说明
[0084]
利用附图对本发明作进一步的说明，但附图中的内容不构成对本发明的任何限
制。
[0085]
图1为本发明的基于物理神经网络的电特性断层成像方法的流程图。
[0086]
图2为全连接神经网络结构图。
[0087]
图3为本发明与helmholtz ept方法在duke人脑肿瘤仿真模型上的重建电容率图像的结果。
[0088]
图4籽本发明与helmholtz ept方法在duke人脑肿瘤仿真模型上的重建电导率图像的结果。
具体实施方式
[0089]
结合以下实施例对本发明的技术方案作进一步说明。
[0090]
实施例1。
[0091]
一种基于物理神经网络的电特性断层成像方法，如图1所示，包括如下步骤：
[0092]
步骤(1)、根据射频场b，得到发射场b；
[0093]
步骤(2)、利用数值方法求解所述步骤(1)得到的发射场b中每点在空间x轴、y轴和z轴上的一阶导数及二阶导数；
[0094]
步骤(3)、构建全连接神经网络，并对麦克斯韦方程进行组合变换得到损失函数；
[0095]
步骤(4)、使用发射场b、步骤(2)得到的一阶导数和二阶导数构成的数据集及步骤(3)得到的损失函数训练步骤(3)得到的全连接神经网络，得到训练后全连接神经网络，所述训练后全连接神经网络以空间任意一点的发射场b以及通过步骤(2)求出在该点的一阶导数为输入，所述训练后全连接神经网络的输出为该点的电容率及电导率，从而得到空间所有点的电容率及电导率。
[0096]
在步骤(1)中，射频场b为射频线圈产生的，且b＝(b
x
,by,bz)，其中b
x
为射频场在空间上x轴的分量，by为射频场在空间上y轴的分量，bz为射频场在空间上z轴的分量，根据射频场b得到发射场b，存在
[0097]
本发明步骤(2)具体为利用数值方法求解发射场b中每点的一阶导数b
x
、一阶导数by、一阶导数bz、二阶导数b
xx
、二阶导数b
xy
、二阶导数b
xz
、二阶导数b
yz
、二阶导数b
yy
和二阶导数b
zz
；
[0098]
一阶导数b
x
由式(7)表示，
[0099][0100]
一阶导数by由式(8)表示，
[0101][0102]
一阶导数bz由式(9)表示,
[0103][0104]
二阶导数b
xx
由式(10)表示；
[0105][0106]
二阶导数b
xy
由式(11)表示，
[0107][0108]
二阶导数b
xz
由式(12)表示，
[0109][0110]
二阶导数b
yz
由式(13)表示，
[0111][0112]
二阶导数b
yy
由式(14)表示，
[0113][0114]
二阶导数b
zz
由式(15)表示；
[0115][0116]
其中(xk,yk,zk)为第k个点的坐标，dx,dy和dz分别为在空间三个方向上离散化后的网格。
[0117]
本发明的步骤(3)中，损失函数的获得方法步骤具体为：
[0118]
步骤(3.1)、对麦克斯韦方程组进行组合变换，得到式(1)
[0119][0120]
其中εc＝ωε-iσ为复介电常数，i为虚数单位，ε为电容率，σ为电导率，ω为核磁共振设备对应的拉莫尔进动频率，μ为磁导率，为nabla算子，为拉普拉斯算子；
[0121]
步骤(3.2)、令bz为0并对步骤(1.1)的式(1)进行分离和重组，得到mr-ept的核心方程式(2)；
[0122][0123]
其中b
x
为发射场b中在空间上x轴的一阶导数,by为发射场b中在空间上y轴的一阶导数,bz为发射场b中在空间上z轴的一阶导数,ε
c,x
为复介电常数在空间上x轴上的一阶导数，ε
c,y
为复介电常数在空间上y轴上的一阶导数，pde为偏微分方程符号；
[0124]
步骤(3.3)、根据式(2)得到损失函数。
[0125]
所述损失函数由式(3)表示，
[0126][0127]
其中由式(4)表示，
[0128][0129]
其中由式(5)表示，
[0130][0131]
其中，α为平衡和的权重且α》0，corr函数为计算电容率梯度和电导率梯度之间的相关性之间的互相关函数，由式(6)表示，
[0132][0133]
其中为x的均值，为y的均值，n为计算点的个数，k为第k个点。
[0134]
需要说明的是，因为使用鸟笼线圈时，则bz的幅值是远小于x轴、y轴上的分量(b
x
和by)，因此是可以忽略不计的。
[0135]
还需要说明的是，本发明的损失函数由两部分构成，第一部分由ept的核心方程组成，目的是使模型的输出的介电常数能够逼近该方程(即为式(4))。
[0136]
其中损失函数涉及到计算电导率和电容率在空间上的一阶导数，根据链式法则，可表示为如下形式，
[0137][0138][0139][0140]
其中f＝wf
ε-if
σ
，则，
[0141][0142][0143][0144][0145]
由于需要对复数求模，而神经网络的运算是实数运算，因此在具体应用中的做法是分别计算pde的实部pdei和虚部pder，然后再代入中，由于因此可以原损失函数的复数运算转换成实数运算。其中关于介电常数对b以及b的梯度的导数可以由神经网络的自动微分方法很方便地求得。
[0146]
损失函数的第二部分参考图像的匹配算法，这是基于以下事实：在人体中电导率和电容率分布具有相似的结构特征，因此电导率和电容率的梯度信息之间具有很高的相关性。因此本发明计算电容率的梯度与电导率的梯度的互相关值，将这个值作为网络的损失函数，使得网络的输出的电导率和电容率在结构上具有较高的相似性。
[0147]
本发明的全连接神经网络分由网络a和网络b两部分组成。网络a构造了f
ε
函数，所述网络b构造了f
σ
函数，所述网络a和所述网络b输入均为空间上任意一点的发射场b以及该点的一阶导数b
x
、一阶导数by和一阶导数bz，输出分别为该点的电容率和电导率。
[0148]
网络a的输入层的激活函数、所述网络a的中间层的激活函数、所述网络b的输入层的激活函数及所述网络b的输入层的中间层的激活函数均为tanh函数，且tanh函数由式(16)表示，
[0149][0150]
网络a的输出层的激活函数、所述网络b的输出层的激活函数均为sigmoid函数，且sigmoid函数由式(17)表示，
[0151][0152]
其中tanh函数将输入映射到(-1,1)之间，sigmoid函数将输入映射到(0,1)之间。
[0153]
当所述网络a和所述网络b输入均为空间上任意一点的发射场b及该点的一阶导数b
x
、一阶导数by和一阶导数bz时，该点的电容率由式(18)表示，该点的电导率由式(19)表示；
[0154]
ε＝(81*oa 1)ecꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
式(18)；
[0155]
σ＝3*obꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
式(19)；
[0156]
其中ε为电导率，σ为电容率，oa为网络a的输出，ob为网络b的输出。
[0157]
需要说明的是，由于人体组织的电导率和电容率的范围已经有实验测得，在外部磁场为1.5t-7t的范围内，人体组织中电容率的范围为ec-82ec，ec为真空电容率，电导率的
范围为0-3s/m，因此在本发明中利用一个简单的线性变换网络a和网络b的输出映射到这些范围中，得到最终的由式(18)得到的电容率和由式(19)得到的电导率。
[0158]
其中，全连接神经网络的映射由式(20)和式(21)表示，
[0159]
ε(xk,yk,zk)
[0160]
＝f
ε
(b(xk,yk,zk),b
x
(xk,yk,zk),by(xk,yk,zk),bz(xk,yk,zk)；w
ε
)
[0161]
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
式(20)；
[0162]
σ(xk,yk,zk)
[0163]
＝f
σ
(b(xk,yk,zk),b
x
(xk,yk,zk),by(xk,yk,zk),bz(xk,yk,zk)；w
σ
)
[0164]
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
式(21)；
[0165]
其中w
ε
和w
σ
均为神经网络的参数，ε(xk,yk,zk)为空间上坐标为(xk,yk,zk)上的电容率值，σ(xk,yk,zk)为空间上坐标为(xk,yk,zk)上的电导率值。
[0166]
本发明的步骤(4)具体为：
[0167]
步骤(4.1)、使用发射场b、步骤(2)得到的一阶导数和二阶导数构成的数据集全部输入步骤(3)得到的全连接神经网络，得到全连接神经网络的输出值，将全连接神经网络的输出值输入至步骤(3)得到的损失函数，令n＝1并进行前向和反向运算得到数据集的当前平均值平均值更新全连接神经网络，进入步骤(4.2)；
[0168]
步骤(4.2)、令n＝n 1，进入步骤(4.3)；
[0169]
步骤(4.3)、使用发射场b、步骤(2)得到的一阶导数和二阶导数构成的数据集全部输入更新后的全连接神经网络，得到全连接神经网络的输出值，进入步骤(4.4)；
[0170]
步骤(4.4)、将全连接神经网络的输出值输入至步骤(3)得到的损失函数，然后进行前向和反向运算得到数据集的平均值判断与迭代阈值的关系，当大于等迭代阈值则进入步骤(4.5)，否则进入步骤(4.6)；
[0171]
步骤(4.5)、将步骤(4.4)得到的平均值更新全连接神经网络，得到更新后的全连接神经网络并返回步骤(4.2)；
[0172]
步骤(4.6)、将当前的更新全连接神经网络为定义为训练后全连接神经网络。
[0173]
数值方法为有限差分方法。迭代阈值为0.001。网络a和所述网络b均由五层全连接层构成，且每层有256个神经元。
[0174]
本发明的全连接神经网络使用python3.6编写，并使用tensorflow深度学习库构建。
[0175]
实施例2。
[0176]
本实例的目的是结合duke人脑肿瘤仿真模型具体介绍如实施例1的一种基于物理神经网络的电特性断层成像方法的实施步骤，包括如下步骤：
[0177]
步骤(1)、根据射频场b，得到发射场b，同时得到参考的电容率和电导率图像数据；
[0178]
步骤(2)、利用数值方法求解所述步骤(1)得到的发射场b中每点在空间x轴、y轴和z轴上的一阶导数及二阶导数；
[0179]
步骤(3)、构建全连接神经网络，并对麦克斯韦方程进行组合变换得到损失函数，本实施例的损失函数具体选取α＝4e-7。
[0180]
步骤(4)、使用发射场b、步骤(2)得到的一阶导数和二阶导数构成的数据集及步骤(3)得到的损失函数训练步骤(3)得到的全连接神经网络，得到训练后全连接神经网络，本实施例训练神经网络具体为使用的优化器为adam，学习率为1e-6；神经网络权重使用正态分布初始化，标准差为1e-4；对每层神经网络的权重使用权重标准化，加速损失函数收敛；每次训练同时输入13864个点的数据计算。当损失函数不再下降时(具体为时)，全连接神经网络训练完成得到训练后全连接神经网络。本发明的训练后全连接神经网络以空间任意一点的发射场b以及通过步骤(2)求出在该点的一阶导数为输入，所述训练后全连接神经网络的输出为该点的电容率及电导率，从而得到空间所有点的电容率及电导率。
[0181]
本实施例将步骤(4)中的电容率以及电导率图像与基于均匀性假设的ept方法(helmholtz ept)分别在视觉上和定量指标上进行比较。图3为两种方法在电容率重建上的结果，图4为两种方法在电导率重建上的结果。从上述两张图可以看出本发明重建出来的电容率与电导率与仿真的电容率和电导率高度一致，且在视觉上明显优于基于均匀性假设的方法。
[0182]
量化指标选择峰值信噪比以及结构相似性来评价本发明与helmholtz ept方法的结果。峰值信噪比与结构相似性的值越高，说明重建出的图像与参考图像越接近，表1是两种方法得出的结果的量化指标。
[0183]
表1、helmholtz ept和本发明得出的结果的量化指标
[0184][0185]
从表1可见，本发明在指标上也能看出在峰值信噪比与结构相似性上的值明显高于helmholtz ept方法，因此本发明在重建精度上明显优于基于均匀性假设的方法。
[0186]
最后应当说明的是，以上实施例仅用以说明本发明的技术方案而非对本发明保护范围的限制，尽管参照较佳实施例对本发明作了详细说明，本领域的普通技术人员应当理解，可以对本发明技术方案进行修改或者等同替换，而不脱离本发明技术方案的实质和范围。