本文旨在提出一種基于頻譜約束的多頻動態電阻抗斷層成像（EIT）算法。在已知成像域內各組分電導率頻譜的情況下，通過重構獨立于頻率的參數——體積分數變化，同時利用多個激勵頻率下的測量電壓差重構一幀時差圖像，從而大大增加測量數據量以改善逆問題的病態性。數值仿真實驗顯示，該算法較傳統阻尼最小二乘算法具有更小的圖像偽影，且在低信噪比情形下具有更小的位置誤差和形變誤差。本研究有望為動態 EIT 提供一種有效利用多頻信息的方法，并為在已知各組分電導率頻譜情況下的動態 EIT 發展提供一個新思路。
引用本文： 曹璐, 楊濱, 李昊庭, 劉學超, 劉本源, 徐燦華, 劉銳崗, 付峰. 一種基于頻譜約束的多頻動態電阻抗斷層成像算法. 生物醫學工程學雜志, 2020, 37(1): 80-86. doi: 10.7507/1001-5515.201901018 復制
引言
電阻抗斷層成像（electrical impedance tomography，EIT）是在人體表面放置一定數量的電極，并通過這些電極向人體施加微小的電流激勵，同時測量響應電壓，使用逆問題求解算法重建出體內阻抗分布及變化信息。它具有無創、無輻射、動態監測、功能成像等優勢，具有良好的應用前景[1]。
EIT 從應用的角度可以分為解剖結構成像和功能成像[2]。解剖結構成像主要有靜態成像和頻差成像。靜態成像是指給定某一時刻邊界電壓值重構體內阻抗分布的絕對值，其不依賴初始阻抗分布但是對模型誤差和測量噪聲敏感；頻差成像則利用兩個頻率下的電壓差值重構體內阻抗隨頻率的變化量，正處于發展初期，具有很強的應用前景。功能成像主要是時差成像，即利用不同時刻邊界電壓差重構體內阻抗隨時間的變化量，時間分辨率高，可以實現連續動態 EIT 監測，而且由于采取邊界電壓差值，很大程度上削弱了模型誤差，是目前在臨床應用最廣泛的方法。但該法存在著依賴初始阻抗分布、對測量噪聲敏感和不適用于阻抗變化大的成像目標等問題[3]。
目前各種 EIT 成像方法面臨的主要問題是病態性，它的實質是物體內部阻抗分布和表面測量電壓之間存在非線性關系，且這種關系難以獲得，導致內部重構阻抗隨著測量電壓中很小的噪聲劇烈變化，非常不利于連續動態的 EIT 監測。實際應用中還需對這種關系進行離散化處理，這使得觀測數據量大大減小，進一步加劇病態性。目前常規解決方法是利用阻抗分布的先驗知識，結合正則化代價函數的方式來得到穩定解[3]。
本文提出一種基于頻譜約束的多頻動態電阻抗斷層成像算法，它同時利用多個頻率下的時差數據重構一幀體積分數變化，成倍增加觀測數據量從而減小逆問題的病態性，且合并正則化方法將提高 EIT 圖像的質量和抗噪能力。該法不同于以往的頻差成像，本質上是一種利用多頻數據的時差成像方式，因此可以實現功能成像。
1 體積分數成像
1.1 體積分數模型
體積分數模型是體內各種組織分布的一種表現形式，通常和有限元法（finite element method，FEM）結合起來，估計體內電導率分布[4]。該模型基于以下三個假設：（1）成像區域是由 T 種具有不同電導率特性的組織構成的，分別是 t1 
tj 
tT。（2）每種組織在不同頻率激勵下的電導率是已知的，j 組織在第 i 個激勵頻率下的電導率是 
。（3）第 n 個剖分面元的電導率等于各組織電導率的線性加權和：
（
）。這其中的加權系數 fnj 即 j 組織在面元 n 中的體積分數。顯然，如果第 n 個面元中僅包含一種組織 j，那它的電導率就是該組織的電導率：σn（ωi）= σij。
以二維區域 Ω 為例，該區域內的電導率分布隨頻率變化定義為 
，其中 x 是空間位置，ω 是激勵電流頻率。然后對 Ω 做有限元剖分，每個面元內的電導率是連續一致的，于是得到了電導率向量 
，其中 N 是剖分面元數。又因為每個面元可能包含 T 種或部分組織，因此在每個面元內對每種組織設一體積分數值，例如 fnj 即面元 n 內第 j 種組織的體積分數，便得到了體積分數矩陣 Fmatrix∈RT*N，其中第 n 列是面元 n 的各組織占比，第 j 行是組織 j 在各個面元的占比，之后將其按列向量化得到 
。最后任意時刻采用 M 個激勵頻率（ω1 
ωi 
ωM）施加激勵電流，每個激勵頻率下可以得到一個邊界電壓向量 
，其中 K 是一個激勵頻率下的所有邊界電壓值個數。將這 M 個邊界電壓向量縱向拼接可得該時刻的邊界電壓向量 υ（t），將該邊界電壓向量和背景幀 υ（t0）作差可以得到該時刻的邊界電壓變化向量?υ（t）。其中 F 形如式（1）：
 |
基于以上假設推導體積分數變化和邊界電壓變化之間的關系即體積分數變化重構的正問題模型。
首先根據上面的假設（3）可以得出體積分數和電導率之間的線性關系：
 |
其中 
是由電導率頻譜構成的系數矩陣，其形式如下：
 |
已知在傳統電導率重構時，不含噪聲理想狀態下的離散正問題：
 |
其中雅克比矩陣 J∈RK*N，它的計算方式有直接求導法和伴隨場法[2]：
 |
結合式（2）和式（4）可以得到理想狀態下體積分數重構的離散正問題方程：
 |
其中的矩陣 J、矩陣 A 和邊界電壓向量 υ 都是隨頻率變化的，
，
，
。只有體積分數向量 F 不隨頻率變化。
1.2 體積分數重構
為了實現同時利用多個頻率下的時差數據重構一幀體積分數變化，構造以下目標函數 Φ，使得所有頻率下的數據差二范數最小。
 |
其中?υ 是任意時刻邊界電壓和背景幀之間的邊界電壓差向量，S 是一個組裝矩陣，內含 M 個元素矩陣，每個元素矩陣是由一個頻率下的雅克比矩陣和系數矩陣相乘得到的，它形如：
 |
該目標函數也可以寫成如式（9）各個頻率下數據差二范數之和的形式，二者是等效的[5]。
 |
之后選擇 Tikhonov 正則化方式[6]，目標函數變成了：
 |
目標函數是可微的，其一階和二階導數由鏈式法則可得。其中 λ 是正則化參數，它代表著正則化程度，維持著數據匹配和噪聲抑制之間的平衡。λ 過大則圖像過正則化，空間分辨率低；λ 過小則圖像正則化不足，主要被噪聲控制，無法分辨重構目標[7]。本文采用 L 曲線法確定正則化參數，該法操作簡單且有效，是目前應用最廣泛的正則化參數選取方法[8]。R 是正則化矩陣，R = I 的形式被稱為 Levenberg-Marquardt 法，簡稱 LM 法；R = diag（STS）的形式稱為標準 Tikhonov 法，也稱 NOSER 法；R 還可以是一些微分算子的離散形式，以起到濾波的作用[9]。本文采用標準 Tikhonov 法。
為了滿足任何一個面元內所有體積分數和為 1，即
，在目標函數中替換 
，對應的 t1 組織通常是背景組織，這樣將重構 T?1 種組織（f2 
fj 
fT）的體積分數變化量?FT?1，而 S 矩陣也將對應地變化成 
。其次為了滿足任何一個體積分數值在 0、1 之間，即 
，對重構結果施加邊界條件以保證 FT?1 中的每個元素被限制在區間［0，1］之間。最終得到目標函數 
：
 '/> |
對?FT?1 的求解將采取逆問題算法，當得到求解式之后，將?υ 代入求解式即可得到?FT?1。
具體流程如下：
步驟 1：設定背景幀 
值。
步驟 2：令一階導數為零得到目標函數的步進方向 
。
步驟 3：利用前一步得到的步進方向求得前景幀的 
值，
= 
+ d，并對其施加限制條件得到 FT?1，具體的限制條件是 
每個元素要 ≥ 0 且 ≤ 1。
步驟 4：利用 FT?1 減去背景幀 
值即可得到區域內部體積分數變化量?FT?1。
最后，為了與傳統時差圖像作對比，利用體積分數和電導率的線性關系將成像區域內體積分數變化轉化得電導率變化，將電導率變化利用色階表示在成像區域內，得到電導率變化圖像。
2 圖像評價指標
參考 Adler 等[10]提出的 EIT 圖像評價體系后，本文采用位置誤差（position error，PE）、形變誤差（shape deformation，SD）、圖像噪聲（image noise，IN）三項指標來評價算法性能。首先定義重建目標區域（region of perturbation，RP），其滿足如下兩個條件[4]：① t2 組織體積分數變化大于最大 t2 組織體積分數變化的 0.5 倍的面元集合；② 按照條件① 提取的所有區域中，覆蓋面積最大的連通區域集合。
2.1 位置誤差
位置誤差衡量重建圖像準確表達真實目標位置的程度：
 |
其中，dRP 和 dREAL 分別表示重建目標重心和預設目標重心至模型中心的距離，l 為重建模型長軸長度。PE 應盡可能的小，并在全局范圍內具有較高一致性，該指標對準確判斷目標的位置起著至關重要的作用。
2.2 形變誤差
形變誤差衡量重建圖像準確表達真實目標形狀的程度：
 |
其中，AC 和 ARP 分別表示區域 C 和區域 RP 所包含的面元個數。SD 應該是小且一致的，較大的 SD 可能會導致對圖像的不正確解釋。
2.3 圖像噪聲
圖像噪聲是衡量算法抗噪聲性能的一項指標：
 |
分子是位于重建目標區域外的體積分數變化的平均值；分母是所有面元體積分數變化的平均值。IN 越小，代表圖像偽影越少。
3 仿真實驗
為了測試基于頻譜約束的多頻動態 EIT 算法的性能，采取圓型 FEM 模型進行仿真實驗。
3.1 實驗步驟
實驗采用 16 電極，對側激勵、鄰近測量模式。如圖 1 所示，生成模板采用 10 層剖分的圓域模型，內含 441 個結點、800 個單元；重構模板采用 8 層剖分的圓域模型，內含 289 個結點、512 個單元。生成和重構采用稀疏不同模板的目的是避免逆問題陷阱[2]。如圖 2 所示，將背景幀設置成均勻分布的 t1 組織，前景幀設置成均勻分布的 t1 組織內含一個由 t2 組織構成的重建目標，前景幀一共設置 5 個，分別具有不同徑向位置的目標 0～4。其中 t1 組織仿真正常腦組織，t2 組織仿真腦缺血組織，電導率頻譜如圖 3 所示，其值通過查閱文獻后截取三個頻點來獲得，截取的方法在討論中詳細闡述[11]。然后用生成模板生成邊界電壓差數據并對其添加兩種水平的高斯噪聲，信噪比分別是 60 dB 和 80 dB。最后分別采用阻尼最小二乘算法和基于頻譜約束的多頻動態 EIT 算法對這 5 個目標進行圖像重建。阻尼最小二乘算法是目前臨床應用最廣泛且最穩定的時差成像算法[12]，它的重構公式是
，其中的 J 是雅克比矩陣[13]，λ 是正則化參數且同樣采用 L 曲線法確定，R 是正則化矩陣取 R = diag（JTJ），最后由一個激勵頻率下的邊界電壓差值?υ 可以算出此頻率下的電導率變化?σ。阻尼最小二乘算法采用頻率標簽 2 下的邊界電壓差數據重構，原因是該頻率下的缺血組織與正常腦組織具有最大的電導率差異。
圖1
有限元剖分示意圖
Figure1.
Finite element meshes
圖2
背景幀和前景幀設置
Figure2.
Background frame and foreground frames
圖3
兩種組織的電導率頻譜
Figure3.
Spectrums of two tissues
本研究還進行了初步的多組織仿真實驗，將組織數 T 擴展到 3 以進一步驗證算法的有效性。其中 t1 組織是背景組織，t2 組織仍是缺血組織，新增一種 t3 組織血液，在此基礎上進行體積分數變化的重構。
3.2 實驗結果
阻尼最小二乘算法和基于頻譜約束的多頻動態算法在信噪比 80 dB 和 60 dB 下的重建結果分別如圖 4 和圖 5 所示。
圖4
信噪比 80 dB 下仿真實驗結果
Figure4.
Reconstructed images under signal-to-noise ratio of 80 dB
圖5
信噪比 60 dB 下仿真實驗結果
Figure5.
Reconstructed images under signal-to-noise ratio of 60 dB
由圖 4 可知，兩種算法在高信噪比的時候，重構目標類似且形狀清晰位置準確，但是傳統的阻尼最小二乘算法仍然在背景部分具有明顯的偽影，尤其是在重構目標靠近中心時。說明當測量情況理想時，兩算法都可以重構出滿意的圖像，但基于頻譜約束的多頻動態算法具有更強的偽影抑制能力。
由圖 5 可知，兩種算法在低信噪比的時候，重構圖像存在明顯不同。傳統阻尼最小二乘算法圖像受噪聲的影響大，圖像偽影嚴重而且目標形變嚴重，尤其是在目標靠近中心時圖像偽影會影響對目標的判斷；而基于阻抗頻譜約束的多頻算法圖像偽影明顯少于前者，且具有更小的圖像形變，即使是在目標靠近中心時仍然可以清晰地判斷目標。說明當測量情況不理想時，阻尼最小二乘算法存在失效風險，而基于阻抗頻譜約束的多頻算法仍然具有良好的成像結果。
綜合圖 4 和圖 5 可知，基于阻抗頻譜約束的多頻算法在不同信噪比下都具有更強的偽影抑制能力，在低信噪比下還具有更小的形變誤差，較阻尼最小二乘算法更適合臨床應用。
算法性能評價結果如表 1～表 3 所示，包括兩種算法在不同位置和不同信噪比下的位置誤差、形變誤差和圖像噪聲。在高信噪比下，基于頻譜約束的多頻動態算法具有同阻尼最小二乘算法相近的位置誤差和形變誤差，且二者的位置誤差都很小，新算法主要優勢是具有更小的圖像噪聲，這點也與重建圖像相對應。在低信噪比下這種優勢更加明顯，基于頻譜約束的多頻動態算法具有更強的噪聲抑制能力，同時具有更小的位置誤差和形變誤差。
表1
位置誤差評價結果
Table1.
Position error evaluation results of two algorithms
表2
形變誤差評價結果
Table2.
Shape deformation evaluation results of two algorithms
表3
圖像噪聲評價結果
Table3.
Image noise evaluation results of two algorithms
多組織的仿真實驗結果如圖 6 所示，該法在成像前就對組織進行了劃分，成像后可以方便地篩選出各個組織的目標圖像，這是傳統阻尼最小二乘算法所沒有的優勢。但缺血目標成像效果不如出血目標，原因是本研究采用的頻點下缺血目標電導率較背景變化遠小于出血目標，因此還需要進一步探討組織個數 T 對成像結果的影響，以實現最優的多組織多目標成像。
圖6
多組織成像目標和重構結果
Figure6.
Multi-tissue case setup and reconstruction results
4 討論
本文在正則化基礎上從算法模型入手建立了一種利用多頻信息進行時差成像的算法。算法基于組織的阻抗譜信息，將體積分數模型引入阻抗成像領域，找到了一個獨立于頻率的成像參數——體積分數變化，同時使用多個頻率下的測量值進行一幀時差成像，能有效增加獨立方程數，改善逆問題病態性，從而提高重構圖像質量。仿真結果表明該算法較阻尼最小二乘算法具有更強的偽影抑制能力，在低信噪比下還具有更小的位置誤差和形變誤差，更適合臨床應用場景。
但需要注意的是：（1）頻率數取決于重構模板和激勵測量模式的選擇。頻率數不是越多越好，只要滿足在逆問題方程中，所有頻率下的測量數大于等于未知的體積分數個數即可。過多的頻率數會增加計算時間，而過少的頻率數將不能體現該算法的優越性，文中選擇頻率數為 3 即可。（2）不同頻點下正常與病變組織電導率差異大小會對譜成像算法的性能造成影響，選擇不同組織的電導率差異較大的頻點，將獲得更好的成像效果。本文首先仿真生成所有頻點下的邊界電壓差數據，然后計算各個頻點下邊界電壓差數據的相關系數，選取正交性最大的三組數據來成像。鑒于本文 EIT 硬件系統在 50 kHz 時具有較高的信噪比，頻點選擇在 50 kHz 及其左右展開，截取的頻譜中心頻率設置在 50 kHz，當其他頻點和 50 kHz 下的邊界電壓差數據具有較高的相關性時，將舍棄這個頻點。（3）本法需要事先獲取完整的組織電導率譜信息，可以通過查閱文獻獲得，也可以通過測量獲得。但以上兩種方式獲得的組織頻譜都存在一定的誤差，會影響成像結果。
綜上所述，本文利用體積分數模型，將多頻信息整合進時差成像中，通過增加重構一幀圖像的信息量，來改善圖像質量。但重構結果很大程度上依賴于引入信息的準確度，因此電導率頻譜信息的獲取和篩選至關重要，最優的頻點選擇方法和頻譜誤差如何影響重構結果以及如何減小這種影響值得進一步探討。
利益沖突聲明：本文全體作者均聲明不存在利益沖突。
引言
電阻抗斷層成像（electrical impedance tomography，EIT）是在人體表面放置一定數量的電極，并通過這些電極向人體施加微小的電流激勵，同時測量響應電壓，使用逆問題求解算法重建出體內阻抗分布及變化信息。它具有無創、無輻射、動態監測、功能成像等優勢，具有良好的應用前景[1]。
EIT 從應用的角度可以分為解剖結構成像和功能成像[2]。解剖結構成像主要有靜態成像和頻差成像。靜態成像是指給定某一時刻邊界電壓值重構體內阻抗分布的絕對值，其不依賴初始阻抗分布但是對模型誤差和測量噪聲敏感；頻差成像則利用兩個頻率下的電壓差值重構體內阻抗隨頻率的變化量，正處于發展初期，具有很強的應用前景。功能成像主要是時差成像，即利用不同時刻邊界電壓差重構體內阻抗隨時間的變化量，時間分辨率高，可以實現連續動態 EIT 監測，而且由于采取邊界電壓差值，很大程度上削弱了模型誤差，是目前在臨床應用最廣泛的方法。但該法存在著依賴初始阻抗分布、對測量噪聲敏感和不適用于阻抗變化大的成像目標等問題[3]。
目前各種 EIT 成像方法面臨的主要問題是病態性，它的實質是物體內部阻抗分布和表面測量電壓之間存在非線性關系，且這種關系難以獲得，導致內部重構阻抗隨著測量電壓中很小的噪聲劇烈變化，非常不利于連續動態的 EIT 監測。實際應用中還需對這種關系進行離散化處理，這使得觀測數據量大大減小，進一步加劇病態性。目前常規解決方法是利用阻抗分布的先驗知識，結合正則化代價函數的方式來得到穩定解[3]。
本文提出一種基于頻譜約束的多頻動態電阻抗斷層成像算法，它同時利用多個頻率下的時差數據重構一幀體積分數變化，成倍增加觀測數據量從而減小逆問題的病態性，且合并正則化方法將提高 EIT 圖像的質量和抗噪能力。該法不同于以往的頻差成像，本質上是一種利用多頻數據的時差成像方式，因此可以實現功能成像。
1 體積分數成像
1.1 體積分數模型
體積分數模型是體內各種組織分布的一種表現形式，通常和有限元法（finite element method，FEM）結合起來，估計體內電導率分布[4]。該模型基于以下三個假設：（1）成像區域是由 T 種具有不同電導率特性的組織構成的，分別是 t1 tj tT。（2）每種組織在不同頻率激勵下的電導率是已知的，j 組織在第 i 個激勵頻率下的電導率是 。（3）第 n 個剖分面元的電導率等于各組織電導率的線性加權和：（）。這其中的加權系數 fnj 即 j 組織在面元 n 中的體積分數。顯然，如果第 n 個面元中僅包含一種組織 j，那它的電導率就是該組織的電導率：σn（ωi）= σij。
以二維區域 Ω 為例，該區域內的電導率分布隨頻率變化定義為 ，其中 x 是空間位置，ω 是激勵電流頻率。然后對 Ω 做有限元剖分，每個面元內的電導率是連續一致的，于是得到了電導率向量 ，其中 N 是剖分面元數。又因為每個面元可能包含 T 種或部分組織，因此在每個面元內對每種組織設一體積分數值，例如 fnj 即面元 n 內第 j 種組織的體積分數，便得到了體積分數矩陣 Fmatrix∈RT*N，其中第 n 列是面元 n 的各組織占比，第 j 行是組織 j 在各個面元的占比，之后將其按列向量化得到 。最后任意時刻采用 M 個激勵頻率（ω1 ωi ωM）施加激勵電流，每個激勵頻率下可以得到一個邊界電壓向量 ，其中 K 是一個激勵頻率下的所有邊界電壓值個數。將這 M 個邊界電壓向量縱向拼接可得該時刻的邊界電壓向量 υ（t），將該邊界電壓向量和背景幀 υ（t0）作差可以得到該時刻的邊界電壓變化向量?υ（t）。其中 F 形如式（1）：
|
基于以上假設推導體積分數變化和邊界電壓變化之間的關系即體積分數變化重構的正問題模型。
首先根據上面的假設（3）可以得出體積分數和電導率之間的線性關系：
|
其中 是由電導率頻譜構成的系數矩陣，其形式如下：
|
已知在傳統電導率重構時，不含噪聲理想狀態下的離散正問題：
|
其中雅克比矩陣 J∈RK*N，它的計算方式有直接求導法和伴隨場法[2]：
|
結合式（2）和式（4）可以得到理想狀態下體積分數重構的離散正問題方程：
|
其中的矩陣 J、矩陣 A 和邊界電壓向量 υ 都是隨頻率變化的，，，。只有體積分數向量 F 不隨頻率變化。
1.2 體積分數重構
為了實現同時利用多個頻率下的時差數據重構一幀體積分數變化，構造以下目標函數 Φ，使得所有頻率下的數據差二范數最小。
|
其中?υ 是任意時刻邊界電壓和背景幀之間的邊界電壓差向量，S 是一個組裝矩陣，內含 M 個元素矩陣，每個元素矩陣是由一個頻率下的雅克比矩陣和系數矩陣相乘得到的，它形如：
|
該目標函數也可以寫成如式（9）各個頻率下數據差二范數之和的形式，二者是等效的[5]。
|
之后選擇 Tikhonov 正則化方式[6]，目標函數變成了：
|
目標函數是可微的，其一階和二階導數由鏈式法則可得。其中 λ 是正則化參數，它代表著正則化程度，維持著數據匹配和噪聲抑制之間的平衡。λ 過大則圖像過正則化，空間分辨率低；λ 過小則圖像正則化不足，主要被噪聲控制，無法分辨重構目標[7]。本文采用 L 曲線法確定正則化參數，該法操作簡單且有效，是目前應用最廣泛的正則化參數選取方法[8]。R 是正則化矩陣，R = I 的形式被稱為 Levenberg-Marquardt 法，簡稱 LM 法；R = diag（STS）的形式稱為標準 Tikhonov 法，也稱 NOSER 法；R 還可以是一些微分算子的離散形式，以起到濾波的作用[9]。本文采用標準 Tikhonov 法。
為了滿足任何一個面元內所有體積分數和為 1，即，在目標函數中替換 ，對應的 t1 組織通常是背景組織，這樣將重構 T?1 種組織（f2 fj fT）的體積分數變化量?FT?1，而 S 矩陣也將對應地變化成 。其次為了滿足任何一個體積分數值在 0、1 之間，即 ，對重構結果施加邊界條件以保證 FT?1 中的每個元素被限制在區間［0，1］之間。最終得到目標函數 ：
| '/> |
對?FT?1 的求解將采取逆問題算法，當得到求解式之后，將?υ 代入求解式即可得到?FT?1。
具體流程如下：
步驟 1：設定背景幀 值。
步驟 2：令一階導數為零得到目標函數的步進方向 。
步驟 3：利用前一步得到的步進方向求得前景幀的 值， = + d，并對其施加限制條件得到 FT?1，具體的限制條件是 每個元素要 ≥ 0 且 ≤ 1。
步驟 4：利用 FT?1 減去背景幀 值即可得到區域內部體積分數變化量?FT?1。
最后，為了與傳統時差圖像作對比，利用體積分數和電導率的線性關系將成像區域內體積分數變化轉化得電導率變化，將電導率變化利用色階表示在成像區域內，得到電導率變化圖像。
2 圖像評價指標
參考 Adler 等[10]提出的 EIT 圖像評價體系后，本文采用位置誤差（position error，PE）、形變誤差（shape deformation，SD）、圖像噪聲（image noise，IN）三項指標來評價算法性能。首先定義重建目標區域（region of perturbation，RP），其滿足如下兩個條件[4]：① t2 組織體積分數變化大于最大 t2 組織體積分數變化的 0.5 倍的面元集合；② 按照條件① 提取的所有區域中，覆蓋面積最大的連通區域集合。
2.1 位置誤差
位置誤差衡量重建圖像準確表達真實目標位置的程度：
|
其中，dRP 和 dREAL 分別表示重建目標重心和預設目標重心至模型中心的距離，l 為重建模型長軸長度。PE 應盡可能的小，并在全局范圍內具有較高一致性，該指標對準確判斷目標的位置起著至關重要的作用。
2.2 形變誤差
形變誤差衡量重建圖像準確表達真實目標形狀的程度：
|
其中，AC 和 ARP 分別表示區域 C 和區域 RP 所包含的面元個數。SD 應該是小且一致的，較大的 SD 可能會導致對圖像的不正確解釋。
2.3 圖像噪聲
圖像噪聲是衡量算法抗噪聲性能的一項指標：
|
分子是位于重建目標區域外的體積分數變化的平均值；分母是所有面元體積分數變化的平均值。IN 越小，代表圖像偽影越少。
3 仿真實驗
為了測試基于頻譜約束的多頻動態 EIT 算法的性能，采取圓型 FEM 模型進行仿真實驗。
3.1 實驗步驟
實驗采用 16 電極，對側激勵、鄰近測量模式。如圖 1 所示，生成模板采用 10 層剖分的圓域模型，內含 441 個結點、800 個單元；重構模板采用 8 層剖分的圓域模型，內含 289 個結點、512 個單元。生成和重構采用稀疏不同模板的目的是避免逆問題陷阱[2]。如圖 2 所示，將背景幀設置成均勻分布的 t1 組織，前景幀設置成均勻分布的 t1 組織內含一個由 t2 組織構成的重建目標，前景幀一共設置 5 個，分別具有不同徑向位置的目標 0～4。其中 t1 組織仿真正常腦組織，t2 組織仿真腦缺血組織，電導率頻譜如圖 3 所示，其值通過查閱文獻后截取三個頻點來獲得，截取的方法在討論中詳細闡述[11]。然后用生成模板生成邊界電壓差數據并對其添加兩種水平的高斯噪聲，信噪比分別是 60 dB 和 80 dB。最后分別采用阻尼最小二乘算法和基于頻譜約束的多頻動態 EIT 算法對這 5 個目標進行圖像重建。阻尼最小二乘算法是目前臨床應用最廣泛且最穩定的時差成像算法[12]，它的重構公式是，其中的 J 是雅克比矩陣[13]，λ 是正則化參數且同樣采用 L 曲線法確定，R 是正則化矩陣取 R = diag（JTJ），最后由一個激勵頻率下的邊界電壓差值?υ 可以算出此頻率下的電導率變化?σ。阻尼最小二乘算法采用頻率標簽 2 下的邊界電壓差數據重構，原因是該頻率下的缺血組織與正常腦組織具有最大的電導率差異。
圖1
有限元剖分示意圖
Figure1.
Finite element meshes
圖2
背景幀和前景幀設置
Figure2.
Background frame and foreground frames
圖3
兩種組織的電導率頻譜
Figure3.
Spectrums of two tissues
本研究還進行了初步的多組織仿真實驗，將組織數 T 擴展到 3 以進一步驗證算法的有效性。其中 t1 組織是背景組織，t2 組織仍是缺血組織，新增一種 t3 組織血液，在此基礎上進行體積分數變化的重構。
3.2 實驗結果
阻尼最小二乘算法和基于頻譜約束的多頻動態算法在信噪比 80 dB 和 60 dB 下的重建結果分別如圖 4 和圖 5 所示。
圖4
信噪比 80 dB 下仿真實驗結果
Figure4.
Reconstructed images under signal-to-noise ratio of 80 dB
圖5
信噪比 60 dB 下仿真實驗結果
Figure5.
Reconstructed images under signal-to-noise ratio of 60 dB
由圖 4 可知，兩種算法在高信噪比的時候，重構目標類似且形狀清晰位置準確，但是傳統的阻尼最小二乘算法仍然在背景部分具有明顯的偽影，尤其是在重構目標靠近中心時。說明當測量情況理想時，兩算法都可以重構出滿意的圖像，但基于頻譜約束的多頻動態算法具有更強的偽影抑制能力。
由圖 5 可知，兩種算法在低信噪比的時候，重構圖像存在明顯不同。傳統阻尼最小二乘算法圖像受噪聲的影響大，圖像偽影嚴重而且目標形變嚴重，尤其是在目標靠近中心時圖像偽影會影響對目標的判斷；而基于阻抗頻譜約束的多頻算法圖像偽影明顯少于前者，且具有更小的圖像形變，即使是在目標靠近中心時仍然可以清晰地判斷目標。說明當測量情況不理想時，阻尼最小二乘算法存在失效風險，而基于阻抗頻譜約束的多頻算法仍然具有良好的成像結果。
綜合圖 4 和圖 5 可知，基于阻抗頻譜約束的多頻算法在不同信噪比下都具有更強的偽影抑制能力，在低信噪比下還具有更小的形變誤差，較阻尼最小二乘算法更適合臨床應用。
算法性能評價結果如表 1～表 3 所示，包括兩種算法在不同位置和不同信噪比下的位置誤差、形變誤差和圖像噪聲。在高信噪比下，基于頻譜約束的多頻動態算法具有同阻尼最小二乘算法相近的位置誤差和形變誤差，且二者的位置誤差都很小，新算法主要優勢是具有更小的圖像噪聲，這點也與重建圖像相對應。在低信噪比下這種優勢更加明顯，基于頻譜約束的多頻動態算法具有更強的噪聲抑制能力，同時具有更小的位置誤差和形變誤差。
表1
位置誤差評價結果
Table1.
Position error evaluation results of two algorithms
表2
形變誤差評價結果
Table2.
Shape deformation evaluation results of two algorithms
表3
圖像噪聲評價結果
Table3.
Image noise evaluation results of two algorithms
多組織的仿真實驗結果如圖 6 所示，該法在成像前就對組織進行了劃分，成像后可以方便地篩選出各個組織的目標圖像，這是傳統阻尼最小二乘算法所沒有的優勢。但缺血目標成像效果不如出血目標，原因是本研究采用的頻點下缺血目標電導率較背景變化遠小于出血目標，因此還需要進一步探討組織個數 T 對成像結果的影響，以實現最優的多組織多目標成像。
圖6
多組織成像目標和重構結果
Figure6.
Multi-tissue case setup and reconstruction results
4 討論
本文在正則化基礎上從算法模型入手建立了一種利用多頻信息進行時差成像的算法。算法基于組織的阻抗譜信息，將體積分數模型引入阻抗成像領域，找到了一個獨立于頻率的成像參數——體積分數變化，同時使用多個頻率下的測量值進行一幀時差成像，能有效增加獨立方程數，改善逆問題病態性，從而提高重構圖像質量。仿真結果表明該算法較阻尼最小二乘算法具有更強的偽影抑制能力，在低信噪比下還具有更小的位置誤差和形變誤差，更適合臨床應用場景。
但需要注意的是：（1）頻率數取決于重構模板和激勵測量模式的選擇。頻率數不是越多越好，只要滿足在逆問題方程中，所有頻率下的測量數大于等于未知的體積分數個數即可。過多的頻率數會增加計算時間，而過少的頻率數將不能體現該算法的優越性，文中選擇頻率數為 3 即可。（2）不同頻點下正常與病變組織電導率差異大小會對譜成像算法的性能造成影響，選擇不同組織的電導率差異較大的頻點，將獲得更好的成像效果。本文首先仿真生成所有頻點下的邊界電壓差數據，然后計算各個頻點下邊界電壓差數據的相關系數，選取正交性最大的三組數據來成像。鑒于本文 EIT 硬件系統在 50 kHz 時具有較高的信噪比，頻點選擇在 50 kHz 及其左右展開，截取的頻譜中心頻率設置在 50 kHz，當其他頻點和 50 kHz 下的邊界電壓差數據具有較高的相關性時，將舍棄這個頻點。（3）本法需要事先獲取完整的組織電導率譜信息，可以通過查閱文獻獲得，也可以通過測量獲得。但以上兩種方式獲得的組織頻譜都存在一定的誤差，會影響成像結果。
綜上所述，本文利用體積分數模型，將多頻信息整合進時差成像中，通過增加重構一幀圖像的信息量，來改善圖像質量。但重構結果很大程度上依賴于引入信息的準確度，因此電導率頻譜信息的獲取和篩選至關重要，最優的頻點選擇方法和頻譜誤差如何影響重構結果以及如何減小這種影響值得進一步探討。
利益沖突聲明：本文全體作者均聲明不存在利益沖突。
