為了探索磁探測電阻抗成像應用于肺部呼吸監測中的可行性，針對新型磁探測電阻抗成像技術的正問題，以真實人體肺部數據構建成像體物理模型，采用了有限元方法得到成像體內部的電勢和電流密度分布，然后根據畢奧-薩伐爾定律分別獲得呼氣末與吸氣末成像體外部的磁感應強度仿真數據。結果表明，吸氣末時成像體周圍的磁感應強度值比呼氣末小 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%，不僅為今后的逆問題圖像重建奠定了基礎，而且對于臨床疾病檢測與判斷也有很好的參考意義。
本文結果為磁探測電阻抗成像應用于肺部呼吸監測仿真研究奠定了基礎，接下來將對電導率重構進行研究。另外，片電極應用于磁探測電阻抗成像還有很多問題有待研究：①以何種形式向成像體注入電流可以獲取更多的數據；②人體的很多組織是各向異性的，頻率會決定組織各項異性的程度，如果考慮此因素，電場強度、電流密度等計算值會有很大差異。因此基于以上原因，將磁探測電阻抗成像應用到肺部呼吸成像，還有很大的研究空間。
引言
隨著科技的進步與發展，電磁場技術在生物醫學研究中得到了廣泛應用，例如對電磁仿生防護神經網絡的研究[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%，不僅為今后的逆問題圖像重建奠定了基礎，而且對于臨床疾病檢測與判斷也有很好的參考意義。
本文結果為磁探測電阻抗成像應用于肺部呼吸監測仿真研究奠定了基礎，接下來將對電導率重構進行研究。另外，片電極應用于磁探測電阻抗成像還有很多問題有待研究：①以何種形式向成像體注入電流可以獲取更多的數據；②人體的很多組織是各向異性的，頻率會決定組織各項異性的程度，如果考慮此因素，電場強度、電流密度等計算值會有很大差異。因此基于以上原因，將磁探測電阻抗成像應用到肺部呼吸成像，還有很大的研究空間。
