2019, Vol. 40
2. 东北大学 中荷生物医学与信息工程学院, 辽宁 沈阳 110169;
3. 东北大学 智能工业数据解析与优化教育部重点实验室, 辽宁 沈阳 110819
2. School of Sino-Dutch Biomedical & Information Engineering, Northeastern University, Shenyang 110169, China;
3. Key Laboratory of Intelligent Industrial Data Analysis and Optimization, Ministry of Education, Northeastern University, Shenyang 110819, China
磁感应断层成像(magnetic induction tomo-graphy, MIT)技术是一种新兴的、非接触的阻抗断层成像技术[1-5].它在测量中不受体表电极的限制, 且系统所产生的磁场容易穿透如颅骨等人体组织, 适合如脑水肿、血肿等疾病的功能诊断.MIT对人体无辐射、成本低、无需植入人体药物, 方便携带且易于长时间监测, 有着重要的研究价值和广泛的医学应用前景[6].
MIT技术根据步骤可分为正问题和逆问题.正问题即已知成像区域内的电导率分布, 通过测量或计算得到检测线圈中测量信号及灵敏度矩阵; 逆问题则根据检测线圈中测量信号变化和灵敏度矩阵, 通过图像重建算法对被测区域电导率分布变化进行成像, 因此这个过程也被称为图像重建过程.现有的图像重建算法主要包括线性反投影(linear back projection, LBP)算法[7]、牛顿-拉夫逊(Newton-Raphson, NR)迭代算法[8]等, 这些算法都是在认为灵敏矩阵合理性的基础上用不同方法进行重建, 改善重建质量.
在图像重建的过程中, 灵敏度矩阵S本身也是影响图像重建质量的重要因素之一.MIT现有灵敏度矩阵存在高度病态问题, 是图像质量提高的障碍.MIT高度病态性表现为灵敏度矩阵对成像区域到检测激励线圈的距离敏感.一方面, 靠近线圈的外部区域与远离线圈的内部区域灵敏度存在很大不同:对边界的电导率高度敏感, 对内部区域的电导率不敏感; 另一方面, 每一种检测激励组合的灵敏度向量, 对检测激励线圈近处的电导率敏感, 对距离远的区域的电导率不敏感.
为了降低灵敏度矩阵的病态性, 颜华等[8]提出一种改进灵敏度矩阵的方法, 用一个专门设计的阈值滤波器增强中心区域的灵敏度, 降低灵敏度矩阵的病态程度.仿真与实验结果表明, 改进的灵敏度矩阵可显著降低正问题的模化误差, 有效提高重建质量, 使重建图像与真实分布更接近.肖理庆等[9]在可有效提高正问题计算精度的有限元模型基础上, 以h细化区域的起始层数、终止层数及三角形有限元内部所插入节点的横坐标、纵坐标为变量, 利用改进粒子群算法h细化有限元模型, 以优化敏感场均匀分布时灵敏度矩阵改善其病态程度.薛倩[10]基于电场中心线的概念提出一种估算灵敏度矩阵的方法, 降低了灵敏度矩阵的病态性, 利用两相流检测实验验证该方法的可行性, 该方法估算的灵敏度矩阵可用于两相流在线监测, 且明显提高敏感场中间区域的空间分辨率.孙青[11]提出将灵敏度矩阵按灵敏区与不灵敏区对灵敏度矩阵进行奇异值修正, 改善灵敏矩阵的病态性, 使灵敏度分布更均匀, 实验证明在一定程度上降低了相对误差.
灵敏度矩阵的病态性本质上受目标物体灵敏度和检测传感器分布影响, 本文主要针对由于检测传感器分布对MIT灵敏度的影响展开, 提出了一种基于聚类优化灵敏度的方法, 将灵敏度向量聚类分组, 利用由能量函数确定的权值进行分组处理, 以减小病态性的影响.与原始灵敏度相比较, 用优化灵敏度能够成功完成重建任务, 且能够有效提高线性反投影等方法的重建效果, 能较好地重建扰动区域, 削弱伪影.
1 MIT理论典型的MIT系统由4部分组成:在目标成像区域周围等间隔排列的线圈阵列、前端电子电路、数据采集系统、计算机主机[12], 如图 1所示.
|
图 1 MIT系统的结构 Fig.1 Structure of MIT system |
磁感应成像问题可以概括为已知检测线圈测量电压数据, 推导生物组织内部电导率, 重建电导率分布图像的过程.MIT的正问题可描述为对涡流问题的求解:
|
(1) |
式中:V为接收线圈的感应电压; σ为所测物体的电导率分布; S为灵敏度矩阵.
标准的图像重建均基于线性近似, 则相应的逆问题形式为
|
(2) |
灵敏度矩阵S通常由正问题求解得到.基于互易定理, 定义灵敏度为
|
(3) |
式中:Ii和Ij分别为线圈i和j作为激励时线圈中的电流; Ai和Aj分别为线圈i和j作为激励的成像区域的矢量磁位; Ωk为有限元单元k的体积; ω为激励电流角频率; m为不同激励-检测测量组合数; n为成像的分辨率.
灵敏度矩阵表示测量电压变化与电导率变化之间的比值关系, 其不同行对应不同激励检测组合的灵敏度向量.灵敏度矩阵为
|
本文针对MIT系统固有结构导致的灵敏度矩阵的病态性, 提出了基于几何差异性的灵敏度矩阵聚类改进方法, 并应用于MIT图像重建算法.
2.1 几何差异性在一种激励检测线圈组合下, 成像中会出现一部分区域由于与线圈组合的距离过远产生不敏感的状况, 这种灵敏度值的区域性异常变小的情况来源于线圈组合的分布, 本文将之称为几何差异性.这种差异性对应于几何的变化, 不能正确反映测量信号与电导率的比, 对根据灵敏度矩阵重建电导率图像是不利的.随着激励检测线圈对应磁感线经过场域的比例变化, 灵敏度向量中不敏感区域占总图像的比例随之改变, 几何差异性会对整个成像区域测量电压变化与电导率变化的对应关系引入误差的大小也作相同变化.而一般灵敏度是将这种不同几何差异性的灵敏度向量一起引入, 直接组成灵敏度矩阵, 几何差异性的引入, 是造成灵敏度矩阵病态性的原因之一, 给重建图像引入对应区域的伪影.
为消除几何差异性对灵敏度矩阵的影响, 本文提出根据不同激励检测线圈组合的敏感区域进行分组.对于一种激励检测组合下, 激励磁场由激励线圈连续发散向四周辐射并激发涡流, 随着距离的增加, 一部分场点的主磁场会急剧减小, 即远离激励线圈的场点的主磁场过小, 进而无法激发足够的次生磁场; 对于已经产生了足够大次生磁场的场点, 其中一部分由于距离检测线圈过远, 使该位置对应次生磁场传播到检测线圈时过小而无法产生足够的检测信号.
基于这两方面的分析, 对于一种检测激励组合与激励线圈和检测线圈的距离都会直接影响该位置次生磁场(对应电压变化值)的检测, 也影响该位置的灵敏度, 本文将两个线圈所夹圆弧(α≤180°)区域判定为靠近两个线圈的区域.在圆形成像区域中, 将两线圈对应圆弧区域设置为敏感区域, 其余成像区域设置为不敏感区域.显然不同的检测激励组合将会得到截然不同的结果, 图 2是8线圈系统的4类组合敏感区域的划分示意图, 粗线包围圆弧区域为敏感区域.它们对应的敏感区域与成像区域的比值分别为1/8, 1/4, 3/8, 1/2, 对应激励线圈与检测线圈的关系依次为邻位、错位、间位、对位.
|
图 2 不同线圈组合下灵敏区域的划分 Fig.2 Division of sensitive areas of different coil combinations (a)—邻位；(b)—错位；(c)—间位；(d)—对位. |
这样处理后, 在n个线圈的MIT系统中, 灵敏度向量总数为n(n-1)/2, 灵敏度矩阵分组数k=n/2.8线圈的灵敏度矩阵分为4组：
|
引入聚类方法中的势能函数, 对MIT系统中根据几何差异性得到的分组进行差异性分析与处理, 消除病态问题.根据不同检测激励组合的灵敏度向量的几何差异, 将灵敏度矩阵分为k组, 并引入组间聚类指标, 分析不同组的几何差异.
2.2 基于几何差异的势能函数计算对于分组组号为c1, c2, c3, c4的数据集, 它的总势能函数为
|
(4) |
式中:E为总势能函数; ei为ci组的势能函数, 反映不同分组的总差异性, 势能越大, 数据的差异性越大.
使用分组势能函数ei评价各个分组的组间差异性.一般情况下分组势能函数的性能评价指标包括内部指标或者外部指标, 这里用外部指标评价组间几何差异性.以扰动分布作为灵敏度分布的参考向量, 选取总偏差量dsum和各组的均值偏差量dcen两个外部指标计算ei:
|
(5) |
|
(6) |
|
(7) |
式中:xij为ci组的第j个灵敏度向量; xi*为ci组的参考向量, 对应扰动分布向量; d(·)为两个向量的欧氏距离; ki为第ci组的向量总数; μi为ci组的均值向量; μi*是ci组参考向量的均值向量.
根据几何差异性的分组结果, 由式(5)可得到每个分组的势能函数ei, 根据势能函数就可以进行几何差异性分析.
对于已知灵敏度矩阵的MIT系统, 本文将根据几何差异性对灵敏度矩阵分组, 对于初始灵敏度为m组向量的MIT系统, 根据几何差异可分为k组, 改进后灵敏度为
|
(8) |
本文提出通过权值处理灵敏度向量分组的方法, 通过权值补偿每一组的灵敏度向量几何差异性造成的不敏感区域的影响.由于势能函数能表征各分组的差异性, 给出根据势能函数计算权值的公式:
|
(9) |
|
(10) |
式中:pi为ci组的权值; ei是ci组的势能函数.
在已知MIT系统中, 基于聚类的优化灵敏度的逆问题求解表示为
|
(11) |
|
(12) |
式中:σ为电导率的分布; σi表示ci组的重建结果; f(·)表示MIT的重建算法; Ci为ci组的灵敏度组合; Vi为ci组对应的测量电压.
基于聚类的优化灵敏度矩阵的MIT重建算法的步骤如下:
1) 建立MIT模型并用有限元分析;
2) 获得灵敏度矩S与归一化的测量电压V;
3) 对灵敏度矩阵进行聚类, 分为k组, 根据式(7)和式(8)获得优化的S*与各组均值向量μi;
4) 根据式(4)~(6)求出各组的势能函数ei;
5) 由式(9)和式(10)求出各组对应的权值pi;
6) 根据式(12)的形式使用重建算法进行重建, 获得重建图像.
3 实验与结果分析 3.1 参数设置为了检验优化灵敏度矩阵的正确性和有效性, 首先将对分组的几何差异性进行验证, 然后实验验证基于聚类优化灵敏度矩阵的重建结果的有效性, 对优化灵敏度前后的重建结果进行对比.实验所使用的电脑处理器为Intel i7-8750H, 内存8 GB, 软件平台为COMSOL 5.3和Matlab 2017b.本文以半径为4.5 cm的圆形作为研究区域, 周围均匀排列着8个线圈.测量线圈的内半径为0.1 cm, 外半径为2.5 cm, 匝数为4, 采用电流激励, 电压检测工作模式, 激励频率为10 MHz.研究区域的背景电导率设置为1 S/m, 病变扰动的电导率设置为2 S/m, 设置了3个圆形扰动, 以区域圆的圆心建立坐标系, 3个圆心坐标分别(2, 2), (-2，2), (-2, -2).
3.2 灵敏度几何差异性实验灵敏度矩阵在不同分组几何差异性存在许多特征, 本节通过分析不同分组的灵敏度向量与扰动的相关性, 以及分组前后的条件数来比较不同分组的几何差异性.
由于灵敏度向量能反应扰动的位置分布, 使用灵敏度向量与扰动分布的相关性来分析分组的几何差异性.为了增强与扰动的对比性, 使用灵敏度向量与无扰动的差值向量计算相关性.本文建立了三扰动模型并进行正问题求解，得到归一化S矩阵, 然后分析不同分组的差值分布与扰动分布的相关性, 确定不同分组之间的灵敏度向量的差异性.通过检验灵敏度矩阵各分组的均值向量分布与向量分布的情况, 并分别与扰动分布进行相关性分析来确定分组之间的几何差异性.
图 3为灵敏度矩阵各分组均值向量的差值分布图.图中蓝色区域的负值表示电导率扰动时该位置的灵敏度相对于无扰动时下降的值, 负值区对应扰动位置.图 4为各分组的负差值分布图.可以看出图中不同组的负差值响应扰动的程度差异显著, 随分组c1~c4的递变, 敏感区域占成像区域的比例依次增加, 灵敏度负差值能越来越好地对应扰动位置.
|
图 3 各分组均值向量的差值分布 Fig.3 Distribution of differences for each group mean vector (a)—c1组；(b)—c2组；(c)—c3组；(d)—c4组. |
|
图 4 各分组的负差值分布 Fig.4 Distribution of negative differences for each group (a)—c1组；(b)—c2组；(c)—c3组；(d)—c4组. |
表 1为各分组的向量分布与扰动分布的相关性.从均值向量的相关性可以得到c2组的相关性仅为不相关, 其他组则都能达到相关.而从分组向量的相关系数可以得到c2组的灵敏度向量是极弱相关, c1组为弱相关, c3组是中度相关, 而c4组则是强相关; 从均值向量的相关性可以得到c2组的相关性仅为不相关, 其他组则都能达到相关.根据几何差异性的分组结果, 分组向量的相关性和图像可以轻易得出不同分组的向量与扰动分布的显著差异性.
|
表 1 各分组的向量分布与扰动分布的相关性 Table 1 Correlation between vector distribution and perturbation distribution of each group |
条件数作为矩阵病态性的关键指标, 在很多病态问题中作为检验指标.图 5对比了3个不同半径的模型分组前后的条件数的变化.由图可以看出优化灵敏度矩阵S*的条件数降低, 这种分组可以降低灵敏度的病态性.
|
图 5 不同模型分组前后的条件数 Fig.5 Number of conditions before and after grouping different models |
由相关性和条件数的分析结果可知, 灵敏度矩阵的几何差异性对灵敏度矩阵响应电导率变化具有显著影响, 基于敏感区域的分组可以很好地反映灵敏度矩阵的几何差异性.
3.3 MIT重建算法验证首先用有限元法求解正问题获得原始灵敏度矩阵, 并分析灵敏度矩阵的几何差异性, 得到聚类优化的灵敏度矩阵, 然后使用LBP和NR迭代两种重建方法进行图像重建, 比较原始灵敏度矩阵和聚类优化灵敏度矩阵在图像重建中的性能, 验证提出方法的有效性.
实验中使用不同的扰动半径进行对比, 半径依次为0.5, 1, 1.5 cm, 使用LBP算法、NR迭代算法进行重建, 并与基于优化灵敏度的线性反投影(pretreatment linear back projection, PLBP)算法和预处理牛顿-拉夫逊(pretreatment Newton-Raphson, PNR)迭代算法进行对比, 重建图像如表 2所示.
|
|
表 2 重建结果 Table 2 Reconstruction images |
采用均方误差和图像相关系数来评价图像重建质量.均方误差是指两个图像对应像素点数值的差异程度, 它被用于表述图像重建的准确性, 公式表示为
|
(13) |
图像相关系数表示两个图像的线性关系, 相关性系数的公式为
|
(14) |
式中:σ为重建的电导率分布向量; σtrue为模型电导率的真实分布向量; σ为重建向量的均值; i表示图像的第i个像素点.
由表 2的重建结果可知, 优化灵敏度矩阵后的重建图像比优化前伪影都减少了.基于优化灵敏度矩阵的LBP算法重建效果得到了改善, 图像的重建效果随着半径增加而变差.基于优化灵敏度矩阵的NR迭代算法重建结果的伪影相对优化前更少.由表 3和表 4可知, LBP算法的均方误差降低了26%~53%, 相关系数提高10%~13%;NR迭代算法的均方误差降低了5%~9%, 相关系数提高4%~8%, 效果良好.
|
|
表 3 不同重建算法的均方误差 Table 3 Mean square error of different algorithms |
|
|
表 4 不同重建算法的图形相关系数 Table 4 Graphic correlation coefficient of different algorithms |
1) 优化灵敏度用于重建效果良好, 能正确反映病变位置与形状, 重建质量有很大提高, 图像误差减小, 伪影减少, 证明了优化方法的有效性.
2) 针对扰动越大, 重建的形变越大, 重建效果越差的情况, 将会从改进重建算法的角度进一步提升重建效果.
由于本文优化法对于成像目标边界锐化程度的提高有限, 而重建结果的高度锐化是医学图像的重要指标之一, 未来MIT另一个研究方向是提高图像的锐化程度.
| [1] |
Muhammad S, Zulkarnay Z, Ibrahim B, et al.
Magnetic induction tomography:a brief review[J]. Journal Teknologi, 2015, 73(3): 91–95.
|
| [2] |
Korjenevsky A, Cherepenin V, Sapetsky S.
Magnetic induction tomography:experimental realization[J]. Physiological Measurement, 2000, 21(1): 89–94.
DOI:10.1088/0967-3334/21/1/311 |
| [3] |
Griffiths H, Stewart W R, Gough W.
Magnetic induction tomography:a measuring system for biological tissues[J]. Annals of the New York Academy of Sciences, 2010, 873(1): 335–345.
|
| [4] |
Soleimani M, Lionheart W R B.
Absolute conductivity reconstruction in magnetic induction tomography using a nonlinear method[J]. IEEE Transactions on Medical Imaging, 2006, 25(12): 15–21.
|
| [5] |
Ma X, Peyton A J, Higson S R, et al.
Hardware and software design for an electromagnetic induction tomography(EMT)system for high contrast metal process applications[J]. Measurement Science and Technology, 2006, 17(1): 111–118.
DOI:10.1088/0957-0233/17/1/018 |
| [6] |
Wei H Y, Soleimani M.
A magnetic induction tomography system for prospective industrial processing applications[J]. Chinese Journal of Chemical Engineering, 2012, 20(2): 406–410.
DOI:10.1016/S1004-9541(12)60404-2 |
| [7] |
Wu X, Yan S, Xu P, et al.
Image reconstruction algorithm for electrical capacitance tomography based on sparsity adaptive compressed sensing[J]. Journal of Electronics & Information Technology, 2018, 14(1): 1250–1256.
|
| [8] |
颜华, 胡丽娟, 王伊凡, 等.
基于改进灵敏度矩阵的ERT图像重建[J]. 仪器仪表学报, 2018, 39(5): 241–243.
( Yan Hua, Hu Li-juan, Wang Yi-fan, et al. ERT image reconstruction based on improved sensitivity matrix[J]. Chinese Journal of Scientific Instrument, 2018, 39(5): 241–243. ) |
| [9] |
肖理庆, 王化祥, 聂文艳.
利用h细化优化电阻层析成像灵敏度矩阵[J]. 天津大学学报(自然科学与工程技术版), 2018, 51(2): 150–158.
( Xiao Li-qing, Wang Hua-xiang, Nie Wen-yan. Optimization of resistance tomography sensitivity matrix using h-refinement[J]. Journal of Tianjin University(Natural Science and Engineering Technology), 2018, 51(2): 150–158. ) |
| [10] |
薛倩.
基于电场中心线的ECT灵敏度矩阵计算[J]. 中南大学学报(自然科学版), 2016, 47(11): 3929–3934.
( Xue Qian. Calculation of ECT sensitivity matrix based on electric field center line[J]. Journal of Central South University(Natural Science), 2016, 47(11): 3929–3934. ) |
| [11] |
孙青.基于电学成像技术的重建算法研究[D].天津: 天津大学, 2011.
( Sun Qing.Research on reconstruction algorithm based on electrical imaging technology[D].Tianjin: Tianjin University, 2011. ) |
| [12] |
Park D, Lee J, Bascunan J, et al.
HTS shim coils energized by a flux pump for the MIT 1.3 GHz LTS/HTS NMR magnet:design, construction, and results of a proof-of-concept prototype[J]. IEEE Transactions on Applied Superconductivity, 2018, 28(3): 1–5.
DOI:10.1109/TASC.2018.2799182 |