为了探索磁探测电阻抗成像应用于肺部呼吸监测中的可行性，针对新型磁探测电阻抗成像技术的正问题，以真实人体肺部数据构建成像体物理模型，采用了有限元方法得到成像体内部的电势和电流密度分布，然后根据毕奥-萨伐尔定律分别获得呼气末与吸气末成像体外部的磁感应强度仿真数据。结果表明，吸气末时成像体周围的磁感应强度值比呼气末小 8.875%。研究结果显示了呼吸时由于肺容积、电导率的不同而导致周围磁感应强度的差异等信息，为后续磁探测电阻抗成像的图像重建以及临床疾病检测奠定了基础。
引用本文： 王慧泉, 尹剑利, 李光旭, 冯彦博, 王金海, 陈瑞娟. 磁探测电阻抗成像呼吸监测仿真研究. 生物医学工程学杂志, 2017, 34(1): 27-33. doi: 10.7507/1001-5515.201604083 复制
引言
随着科技的进步与发展，电磁场技术在生物医学研究中得到了广泛应用，例如对电磁仿生防护神经网络的研究[1]、电磁辐射在物理人建模中的应用[2]、颅内电磁应力分析等[3]。尤其是电子计算机断层扫描（computed tomography，CT）、核磁共振、超声以及核医学等医学成像技术已成为相关疾病的主要检查手段。
人体肺部常见疾病包括慢性支气管炎、哮喘、支气管扩张症、肺炎、肺结核及肺部肿瘤等。用于肺部疾病检查的临床技术有多种，但现有的常规影像学诊断还不是连续成像，难以反映具有特殊周期性、时域性的呼吸过程，而且它们大多只能在出现结构性病变之后才能给出检查结果。基于电特性检测的电阻抗成像技术（electrical impedance tomo-graphy，EIT）则可发现器质性病变之前出现的功能性异常，用于肺部呼吸疾病的早期筛查，现已在医学成像领域取得了快速发展。
1978 年 Henderson 等[4]基于成像原理得到了首幅肺部层析图像，从此肺部呼吸电阻抗成像技术的研究一直都是业内的热点。许多科研机构和商业公司投入了大量的资金进行相应系统和产品的研究和临床应用，并取得了阶段性研究成果。德国德尔格公司（Drager）研发的肺部呼吸机于 2010 年通过国际质量体系认证，开启了电阻抗断层成像临床应用市场化的先河[5]。Ferrario 等[6]利用德尔格公司的肺部无创检测系统（Pulmovista 500）将图像分为 4 个兴趣区，根据呼吸和心跳频率的不同，进行肺部与心脏区域的成像分析。Zhao 等[7]对图像数据采用不同的变量系数进行筛选分析，可以有效提取电阻抗临床数据指标。巴西的迪斯/飞利浦（Dixtal/Phillips）、瑞士的斯威堂穆（Swisstom）、美国的康尔福盛（Carefusion）等公司也投入了大量的资金进行这方面的研究，并申请了若干相关发明专利。在肺部人工通气[8]、肺血管阻塞症[9]、肺灌注显像等领域[10]，肺部电阻抗成像都取得了很好的研究成果。
此外，对于肺部呼吸功能的电阻抗成像，电极激励、成像物理模型、重建图像算法等仍是该领域的研究热点与难点。天津大学、第四军医大学以及河北工业大学等也在激励源硬件系统的设计以及图像重建算法等方面进行了相关实验研究[11-13]。尽管各种电阻抗功能成像方式在应用于肺部呼吸的检测方面都有杰出的研究成果，但是很多都是接触式成像方式，不仅获得的测量信息量小，并且测量精度还受到许多干扰因素的影响，例如：皮肤表面湿度与粗糙程度以及测量时所施加的压力等，对于阻抗频谱以及图像分辨率有很大的影响。
磁探测电阻抗成像（magnetic detection electri-cal impedance tomography，MDEIT）是利用测量线圈对成像体外部一定距离的磁场值进行测量，进而根据相关的成像算法重建成像体内部的电导率分布图像。作为一种不完全接触式成像方式，磁探测电阻抗成像能够减少电极接触等造成的不利因素影响。本文的肺部呼吸过程中周围磁场变化研究就是基于此展开的，并发现由于肺容积、电导率的不同而导致周围磁感应强度的差异等信息，为后续磁探测电阻抗成像的图像重建以及临床疾病检测奠定了基础。
1 基础与关键技术
1.1 磁探测电阻抗成像
Ahlfors 等[14]首次提出利用测量装置检测成像模型外的磁场数据，进而根据场域外的磁场测量值重构成像体内部的电导率分布图像的技术，并将其命名为磁感应成像（magnetic impedance tomogra-phy，MIT）。1999 年，Tozer 等[15]通过实验证明了由成像体外部磁场数据进行二维电场强度、电流密度图像重建的可行性。为了将新技术与磁感应电阻抗成像区分开来，Ireland 等[16]在 2004 年又重新对此技术命名为磁探测电阻抗成像。
磁探测电阻抗成像技术是利用贴在成像体表面的电极向其注入一定频率的交变电流 I，在成像模型内部产生电压 V、电场强度 E 和电流密度 J，进而在成像体场域外部产生感应磁场 B，用磁场检测装置，测量成像体内部电势电流在其周围一定距离产生的磁场，然后基于重建图像的逆问题求解获得产生磁场的电流与电场强度分布，进而重构出成像目标体的电导率图像。磁探测电阻抗成像原理流程图如图 1 所示。
图1
磁探测电阻抗成像原理流程图
Figure1.
Schematic diagram of Magnetic Detection Electrical Impedance Tomography
1.2 正问题描述
将电极贴在成像区域 Ω 的边界∂Ω上，通过电极向成像体注入电流 I，由于电磁感应，在成像体的内部会产生相应的电势分布和电流密度分布。依据麦克斯韦方程以及欧姆定律，拉普拉斯方程适用于区域 Ω 内部的电势分布 φ 和电导率分布 σ，即
| $ \nabla \cdot \sigma \nabla \varphi = 0 $ |
在边界∂Ω上满足纽曼边界条件，有
| $ - \sigma {\rm{ }}\partial \varphi /\partial {n} = {n} \cdot {j} $ |
式中 n 为边界∂Ω上法线单位向量。将式（2）经过一定变换可以得到
| $ - \partial \varphi /\partial {n} = \left\{ \begin{array}{l} I/A + \;\;\;\;{\text{电流流入点}}\\ {\rm{ - }}I/A-\;{\text{电流流出点}}\\ {\rm{0}}\;\;\;\;\;\;\;\;\;\;\;\;\;{\text{边界其他点}} \end{array} \right. $ |
式中A+为电流注入电极的区域，A—为电流流出电极的区域。根据电导率分布 σ 和纽曼条件，求解拉普拉斯方程可以获取电势 φ 的大小，进而成像区域 Ω 内部的电流密度为
| $ {J} = - \sigma \cdot \nabla \varphi $ |
在此基础上，通过毕奥-萨伐尔定律可以求出目标体外部的磁感应强度
| $ B(r') = \frac{{{\text{μ}_0}}}{{4\pi }}\text{∫}_\Omega {{\rm{J(}}r{\rm{)}}} \times \frac{{r' - r}}{{\left| {r' - r} \right|}}{\rm{d}}v $ |
式中：μ0为真空磁导率；r′=（x′，y′，z′）为场点；r=（x，y，z）为源点。
2 模型仿真
2.1 立方体模型
正问题解析解在电导率一定且匀称分布时较易求解，因此首先对电导率均匀分布的立方体成像模型进行了图像仿真实验，进而验证有限元的计算精度。仿真模型如图 2 所示，尺寸为 10 m×10 m×10 m 的立方体。在该模型两侧表面以一对片电极的形式注入电流，频率为 31.25 kHz，对片电极覆盖模型左右两个侧面，大小为 1 A，电流方向沿 X 轴方向正方向。采用有限元方法，将模型剖分成大小为 1 m×1 m×1 m 的小单元模块，总共 10×10×10 个小立方体单元，11×11×11 个节点。假设在模型内部电导率均匀分布为 1 S/m，并以立方体的中心处为零参考电势，计算成像体内的电势分布。
图2
立方体模型及其中间层电势仿真分布图
Figure2.
Cube model and its middle layer electric potential distribution simulation
利用有限元方法分析获得的中间层（第 5 层）的电势分布图像如图 2 所示。该成像体的电势解析解可由下式获得：
| $ \varphi \left( {x,y,z} \right) = - 0.05 + 0.01 \times y $ |
比较成像体内各节点处电势的数值解与解析解的差异，得到其平均误差为 3.589 2×10–6%（即把每个节点的数值解减去式（6）计算出的解析解，差值相加再除以总节点数）。
由式（4）计算求得电流密度，进而通过毕奥-萨伐尔定律可以得到目标成像体周围的磁感应强度分布。
2.2 圆柱体模型
为了分析磁探测电阻抗成像在呼吸过程监控中的可行性，将模型进一步向真实模型靠拢，构建立方体圆柱以模拟整个胸腔体，并进行仿真分析。在 0.5 m×0.5 m×0.5 m 的立方体内部构建半径为 0.25 m 的圆柱模型，如图 3 所示。首先利用有限元方法将立方体剖分成 50×50×10 个单元，再从中选取在圆柱体内部的单元而剔除外部的单元和节点，剩下的即为成像体研究对象（共 10 层 19 760 个单元 22 847 个节点）。在圆柱体模型的两侧通过一对片电极注入电流，片电极为 0.785 m×0.01 m 的长方形薄片，厚度忽略不计。电流频率为 31.25 kHz，强度为 1 A，方向为沿 X 轴正方向。我们假设该模型内部的电导率以 1 S/m 均匀分布。采用与立方体模型一样的方法进行求解，可得到第 5 层单元节点的电势分布图，如图 3 所示，然后即可利用毕奥-萨伐尔定律求得第 5 层单元节点外部的磁感应强度值。
图3
圆柱体模型及其肺部模型的电势分布图
Figure3.
Electric potential distribution of the cylinder model and the lung model
为使成像体模型更加逼真，在此模型的基础上，将模型大小设定为半径为 0.25 m，高为 0.5 m 的有两个小圆柱异质体的圆柱体模拟胸腔肺部模型，即分别在该圆柱体模型内部对称放置两个半径为 0.05 m 的等高圆柱，如图 3 所示，内部两个圆柱电导率为 1 S/m，其余区域电导率以 0.548 S/m 均匀分布。在圆柱体的侧面加上一对片电极，大小为 0.10 m×0.01 m，厚度忽略不计，肺部模型第 5 层单元节点的电势分布图如图 3 所示，并利用上述方法求得第 5 层单元节点外部的磁场分布值。利用该模型进行的仿真实验结果表明，成像体有异质体时，周围的磁感应强度比无异质体均匀分布时的对应值平均大 2.334%。
3 磁探测电阻抗成像在肺部呼吸中的应用
3.1 呼吸过程概述
呼吸是人体重要的生命体征之一。呼吸具有其特殊的周期性与时域性，呼吸全过程示意如图 4 所示，在呼气与吸气时，肺部的气压与体积是有很大的区别的，肺容积在呼吸过程中的呈曲线变化，b 点是吸气时，由于呼吸作用使外界空气进入肺部，肺内气压减少，肺容积增大。人体呼气末，肺容积为 1 500 mL；人体吸气末，肺容积为 5 000 mL。
图4
呼吸过程及肺容积在呼吸过程中的变化曲线
Figure4.
Breathing process sketch in the process of breathing and lung volume change curve
生物体进行新陈代谢时，参与新陈代谢的组织、器官的体积大小会随时间而变化，进而组织阻抗值也会随之改变，由生物组织的血液阻抗可以计算测量出血液速度[17]，进而发现血液阻抗值发生了变化。已有的生理数据表明[18]，人的肺部阻抗在呼气末时为 727 Ω，而在吸气末时为 2 623 Ω，肺部电导率则在 0.042～0.138 S/m 范围内变化。我们根据在 20 kHz～100 kHz 频率下测试的结果，人体组织电阻率参数的平均值和 95% 置信区间值如表 1 所示。
表1
人体各组织器官的电阻率参数表
Table1.
Table of electrical resistivity parameters of human tissue and organs
3.2 肺部真实模型
在以上仿体研究基础上，通过肺癌图像数据库联盟（the lung image database consortium，LIDC）[19]的人体肺部CT图像进一步构建正常人体真实的肺部模型进行仿体实验，成像区域半径为 0.25 m，高 0.5 m，为了仿真数据更为精确，将目标体剖分为 64×64×10 个单元，在目标体的周围贴放一对长方形电极，忽略电极厚度，电极大小为 1 cm×10 cm，电流频率为 31.25 kHz，电流强度大小为 50 mA，方向沿 X 轴正方向，依据上述的正问题的求解方法流程，先对成像体采用有限元法求解电场强度，进而得到电流密度，然后根据毕奥-萨伐尔定律计算得出场域外部的磁感应强度值。本文仿真实验与上述其它仿真实验均是在计算机平台（配置为：CPU 类型：Skylake 双核 i7 处理器 6Y30；显卡类型：集成显卡；内存：4 GB；硬盘容量：500 G）上，借助了 Matlab 软件平台完成的。
在仿体模型求解时，为突出呼吸过程中周围磁场的变化，假设心脏、肝脏和胃这几个异质体的电导率为各向同性且在呼气与吸气时为恒定值。根据文献[20]，具体参数分别设置为 0.126 S/m、0.236 5 S/m 和 0.357 S/m，仿体其余部分设定为肌肉，电导率为 0.584 8 S/m，在有限元计算中，目标体的剖分单元与 CT 扫描体素的大小一致。
通过线圈测量成像目标体周围的磁感应强度，测量平面平行于 XY 平面且位于肺部中央，为便于观察，每个线圈上有 36 个等角度间隔的测量点，每个测量点上的切线方向即为线圈方向。
人体在呼吸时，肺的体积有很大差异，根据 LIDC 数据库中的肺部 CT 图，利用各区域异质体像素的差异进行区分，并构建人吸气时和呼气时的肺部仿真模型如图 5 所示。在呼气末时，由上述分析的数据可知，此时肺容积为 1 500 mL，电导率为 0.137 S/m，按照求解正问题的计算方法，求出仿体模型中间层的电势分布，然后再根据电势求出电场强度 E，进而求出电流密度 J，根据毕奥-萨伐尔定律计算出场点的磁感应强度，进而求出成像体周围磁感应强度分布。
图5
真实模型及其对应的磁感应强度数值解
Figure5.
True model and its corresponding numerical solution of magnetic induction intensity
在吸气末时，肺容积 5 000 mL，电导率为 0.04 S/m，按照与呼气时进行同样的方式进行求解，可以分别求解出相应的电势、电场强度以及周围的磁感应强度分布，对应的分布图像如图 5 所示。将仿体呼气与吸气周围磁感应强度大小进行比较分析，以（0.15，0.61，0.45）、（0.40，0.63，0.45）、（0.15，0.09，0.45）以及（0.40，0.07，0.45）这 4 个场点为例，如表 2 所示。对比呼气末与吸气末两过程中所有场点的测量值，磁感应强度对比如图 6 所示。结果表明，仿体吸气末周围的磁感应强度值比呼气末时的对应值小 8.875%。
表2
4 个场点呼气末和吸气末时周围磁感应强度
Table2.
Magnetic induction intensity of four points outside the object at the ends of exhaling and inhaling
图6
呼气与吸气末两种情况下对应场点的磁感应强度对比
Figure6.
Comparison of the magnetic induction intensity between the two ends of inhalation and exhalation cases at the corresponding sites
4 结语
本文针对磁探测电阻抗成像的正问题，通过均匀立方体、非均匀圆柱体以及模拟肺部模型，采用有限元方法获得电势、电场强度和电流密度的数值解，根据毕奥-萨伐尔定律分别获得呼气末与吸气末时目标体外部的磁感应数据。对均匀立方体模型、非均匀圆柱体模型以及模拟肺部模型进行了仿真实验，结果表明，吸气末时，仿体模型周围的磁感应强度值比呼气末时减小 8.875%，不仅为今后的逆问题图像重建奠定了基础，而且对于临床疾病检测与判断也有很好的参考意义。
本文结果为磁探测电阻抗成像应用于肺部呼吸监测仿真研究奠定了基础，接下来将对电导率重构进行研究。另外，片电极应用于磁探测电阻抗成像还有很多问题有待研究：①以何种形式向成像体注入电流可以获取更多的数据；②人体的很多组织是各向异性的，频率会决定组织各项异性的程度，如果考虑此因素，电场强度、电流密度等计算值会有很大差异。因此基于以上原因，将磁探测电阻抗成像应用到肺部呼吸成像，还有很大的研究空间。
