基于改進牛頓-拉夫遜算法的腦出血磁感應斷層成像研究

曹弘貴 葉 波 姜 瑛 羅思琦 曹眾楷 歐陽俊林

(昆明理工大學信息工程與自動化學院 昆明 650500)(云南省人工智能重點實驗室 昆明 650500)

1 引言

腦出血 (IntraCerebral Hemorrhage, ICH)是由于腦內血管破裂或滲透引發的出血性病變，占腦卒中的15%～20%，其死亡率高達50%，是一種嚴重危及生命的疾病[1]，對于腦出血疾病，早期、準確的檢測及診斷對挽救患者生命至關重要，后面持續開展對腦出血病變的監測對于有效治療也具有重大意義[2]。

傳統的醫學影像檢測技術有計算機斷層成像(Computed Tomography, CT)、磁共振成像(Magnetic Resonance Imaging, MRI)和光學相干層析成像 (Optical Coherence Tomography,OCT)[3,4]。CT重建圖像的分辨率較高但病灶和背景間對比度低，且輻射劑量較大，MRI能對軟組織成像且無放射性損害，但偽影也較CT多[5]，OCT多用于眼底疾病的成像檢測，也被證明其用于腦出血檢測的可行性，但其穿透深度較低無法檢測深層的顱腦疾病[6]。此外，CT和MRI設備龐大，成本高、成像時間長，在緊急情況下應用受到限制，無法給患者提供快速的檢測和診斷及實時和持續的監測[3]。磁感應斷層成像 (Magnetic Induction Tomography, MIT)是一種基于電磁檢測原理對人體的被動電磁特性進行成像的功能性成像技術，具有非接觸、安全無創、低成本、便攜、無需耦合劑、可持續監測等特點[7]。MIT因其產生的激勵磁場能夠穿透低電導率的頭皮和顱骨，在醫學成像領域有著巨大的潛在優勢，能實現對腦出血的早期診斷和床旁實時、長期、無創的監測[8]。因此，將MIT用于腦出血成像對腦出血的治療、監測具有重大意義。

近年來，很多學者對腦出血MIT開展了研究。由于顱腦組織的復雜結構及其各向異性使正問題的求解變得十分困難，目前在腦出血MIT建模與仿真計算的研究中多使用簡化的腦模型，如2維同心圓、3維同心球及圓柱模型等[9-12]，與實際頭部幾何結構相差較大，在求解正問題時不可避免會帶來相應的計算誤差，導致后續圖像重建結果誤差增大；在MIT圖像重建的研究中，柯麗等人[13,14]建立包含頭皮、顱骨、腦實質、腦出血的顱腦模型，并在其基礎上進行了腦組織的阻抗分析，利用改進反投影算法進行圖像重建，文獻[8,15,16]對多頻成像方法及線圈陣列方式進行了研究，并采用Tikhonov算法進行圖像重建，Zhang等人[17]提出了一種自適應閾值的分裂布雷格曼 (Split Bregman, SB)算法，降低了重建圖像的噪聲和面積誤差，但并未分析重建圖像的相關系數且最低只能對4 ml的腦出血進行成像，Chen等人[18]提出了一種自適應閾值優化的Tikhonov算法，實現了不同出血量腦出血的成像，但重建圖像的相關系數有待提高，周曦等人[19]從初始值和正則化因子選取、靈敏度矩陣更新策略以及圖像閾值優化等4個方面對牛頓-拉夫遜 (Newton Raphson, NR)算法進行改進，實現了不同噪聲的圖像重建，但重建圖像質量還有待提高，Han等人[12]在NR算法目標函數中加入加權矩陣和L1范數懲罰項提高了對腦出血圖像重建的精度，但增大了計算量，算法收斂較慢。從上述研究可以看出，線性反投影 (Linear Back Projection, LBP)算法和Tikhonov正則化算法等單步圖像重建算法可直接求解場域內的電磁特性分布，成像速度較快，但圖像重建質量相對較低；Landweber算法、NR算法和共軛梯度最小二乘法 (Conjugate Gradients Least Squares, CGLS)等迭代圖像重建算法需不斷迭代求解場域內的電磁特性分布，成像速度比單步成像算法慢，但圖像重建質量相對較高，在迭代圖像重建算法中，NR算法應用最為普遍，其成像質量相對較高，收斂性相對其他迭代圖像重建算法較好，但在腦出血MIT中由于組織的電導率極小，發生腦出血時測量數據較小的變化會對重建結果產生較大的影響，導致算法的穩定性和收斂性較差，目標與背景間偽影較大，圖像重建質量也有待提高，因此有必要對其進行改進，提高其穩定性、收斂性和重建質量，同時高質量的重建圖像也為后續的醫學圖像處理及臨床中進行輔助診斷和評估患者預后提供基礎[20]。

針對上述問題，為提高解的穩定性和重建圖像的質量，本文提出一種改進的NR算法，在目標函數中加入自適應加速正則化懲罰項和L2范數正則化懲罰項，并引入投影算子P施加物理意義上的約束，提高算法每一步迭代的效率，減少重建圖像的偽影，改善成像質量；其次為驗證算法在腦出血MIT中的實用性以及提高腦模型的真實性，構建了包含頭皮、顱骨、腦脊液和腦實質4種組織的真實3維顱腦模型，并開展了3個位置處出血量各為24 ml,14 ml, 2 ml的腦出血仿真實驗，通過模型計算得到圖像重建所需的相位差檢測值和靈敏度矩陣，然后利用改進NR算法與Tikhonov, Landweber, NR,CGLS、優化迭代策略的NR等算法分別對設置的腦出血進行磁感應斷層成像，驗證所提算法的有效性。

2 腦出血MIT檢測原理

腦出血MIT檢測的基本原理：采用電磁感應原理，向激勵線圈中通入正弦交流電，激勵線圈在被測區域內產生初級交變磁場B0，初級交變磁場B0作用于人腦時，會在腦組織中感應出渦流，從而產生次級磁場 ΔB，在不同病理狀態下次級磁場ΔB的值會發生改變，可通過檢測線圈得到初級磁場B0和疊加磁場B0+ΔB之間的相位差 Δφ，再通過圖像重建算法得到成像域內的電導率分布圖像，進而通過電導率分布圖像診斷腦出血的嚴重程度。MIT的基本原理及初、次磁場關系如圖1所示。

圖1 MIT的基本原理及初、次磁場間關系

對于較小的電導率擾動 Δσ，Δφ和 Δσ間的關系可以簡化為近似離散線性方程[21]

其中，ΔφM×1∈R是 線圈測量相位差，SM×N ∈R是靈敏度矩陣或雅可比矩陣，表示成像區域內電導率分布變化引起的線圈檢測相位差變化率[21]，ΔσN×1∈R是 被測區域內的電導率分布變化，M是激勵-檢測線圈組合數量，N是被測域內的網格數量，v是噪聲。

腦出血MIT是將通過模型或實驗獲得的檢測值Δφ和成像域內的靈敏度矩陣S求解未知的電導率特性分布 Δσ，但是由于MIT中線圈數量有限，因此線圈檢測值 Δφ的數量遠遠小于待求解的電導率特性分布 Δσ的數量，同時 Δφ的微小變化將對Δσ的值產生較大影響，因此式(1)是一個欠定、病態的方程，需要通過一定的圖像重建算法進行求解，將求解得到的 Δσ按幅值大小可繪制成反映腦出血情況的圖像。

3 改進的Newton-Raphson圖像重建算法

腦出血MIT圖像重建是根據檢測信號和場域內的靈敏度分布重建頭部不同組織中的電導率特性分布。NR算法的基本思想是不斷迭代更新電導率分布值，尋找電導率分布的最優解，使得目標函數的值最小[22]。采用NR算法求解MIT圖像重建問題時通常以相位差測量值和計算值的誤差范數平方作目標函數

其中，F(Δσ) 和 Δφ分別是相位差的計算值和測量值，NR算法的迭代公式為

其中，Δσk+1是待求解的成像區域內電導率分布變化，?f(Δσk)是 梯度方向，Hk是Hessian矩陣。

NR算法每一步迭代都需要求解目標函數的2階偏導數，計算復雜度較大，同時Hesse矩陣具有較嚴重的病態性，尤其在腦出血MIT時生物組織的電導率極小，腦出血測量數據較小的變化會對重建結果產生較大的影響，導致算法的穩定性和收斂性較差，重建圖像中背景與腦出血間的偽影較大。

采用優化方法求解式(2)的過程是不穩定的，通過在目標函數中增加L2范數正則化懲罰項來降低解的不穩定性，同時為了提高NR算法中每步迭代的計算效率，提高圖像重建質量，在目標函數中加入與前一刻計算結果相關的自適應加速正則化懲罰項進行約束，加快算法的收斂速度，該懲罰項在目標函數中的權重隨著迭代次數的增加而自適應改變，改進算法的目標函數為

其中，F(Δσ) 是相位差計算值，Δφ是相位差檢測值，βk是自適應加權系數，隨著迭代次數的增加而自適應改變該懲罰項的權重，Δσk是當前時刻電導率分布變化，Δσk-1是前一時刻電導率分布變化，α ＞0，是用于控制解平滑性和連續性的正則化參數，從理論上確定其值較為困難，一般根據經驗進行設定，R是正則化矩陣，一般為單位矩陣。將式(4)展開并求其下降梯度?f(Δσk)和 Hesse矩陣Hk

其中，F′(Δσ)是 靈敏度矩陣S，將F′(Δσ)用S表示，并將式(5)、式(6)代入式(3)后得到

其中，βk除了可以根據經驗選取為固定值，還可以利用v方法求解其在每次迭代中的最優值，v方法是由Brakhage提出的一種用于求解不適定問題或條件數很差的適定問題的計算方法，與其他用于求解此類問題的方法相比，采用v方法計算目標函數誤差項的加權系數，僅需較少的計算步驟，便可得到較小的最小誤差范數[23]。計算公式為

為了進一步加快算法的收斂，將LBP算法計算得到的電導率分布值作為改進算法的迭代初值，同時根據腦出血MIT檢測的實際情況在改進算法中引入投影算子P確保解收斂。綜上，所提改進NR算法的迭代公式如式(9)，其流程圖如圖2所示。

圖2 改進NR算法流程圖

其中，LBP算法的計算結果作為改進NR算法的迭代初值；P是非負凸集上的投影算子，通過投影算子P引入物理意義上的約束，確保每次迭代的解都非負有界且收斂于一個凸集上；βk使用v方法確定其最優值，v的值一般通過經驗法選取[23]，本文v取8。在迭代過程中設置迭代終止條件為|Δσk+1-Δσk|＜ε或k ＞kmax，當迭代過程中兩個相鄰解間的差小于設定的較小閾值ε=1×10-12時停止迭代，同時設置了最大迭代次數kmax=1 000，防止算法陷入半收斂或無限迭代，進行重復的無效計算，占用計算資源。滿足迭代終止條件時輸出電導率分布Δσk+1，并按其幅值大小繪制得到重建圖像。

4 建模與仿真

4.1 腦出血MIT有限元模型的建立

構建具有真實幾何形狀的顱腦模型，能提高對模型渦流場、磁場的求解精度，在此基礎上開展線圈檢測值和場域內靈敏度的仿真計算，為后續的圖像重建提供可靠的仿真實驗結果。建立的仿真模型共包含3個域：真實顱腦模型、空氣域、線圈域。顱腦模型的幾何結構采用電氣與電子工程師協會(Institute of Electrical and Electronics Engineers,IEEE)、國際電工委員會 (International Electrotechnical Commission, IEC)和歐洲電工標準化委員會 (European Committee for Electrotechnical Standardization, CENELEC)提供的標準比吸收率(Specific Absorption Rate, SAR)測量規范中相同的幾何結構，利用COMSOL Multiphysics對MRI數據進行放樣、轉換、旋轉、縮放等操作后得到3維頭部輪廓；將具有類似介電特性的組織進行合并，保留頭皮、顱骨、腦脊液和腦實質等4種主要頭部組織[24]，并將得到的3維頭部輪廓按各組織1:0.92:0.87:0.84的比例進行縮放[25]，最終得到具有4種組織的真實3維顱腦模型，空間中笛卡爾坐標系的原點位于建立的腦模型的幾何中心。

為模擬實際檢測環境，設置一個足夠大的球型空氣域，半徑為70 cm；根據線圈檢測特性和顱腦模型的結構，采用8線圈均勻布置的方式，依次編號1～8，每個線圈大小、形狀完全一致，均勻分布于腦模型中與頭皮距離為12 cm的圓上，線圈外徑為20 mm，內徑為15 mm，高為30 mm，匝數為100，線圈導線電導率為 6×107S/m，導線的截線面積定義為 10-7m2，建立的腦出血MIT有限元模型如圖3所示。

圖3 腦出血MIT有限元模型

將模型中線圈的材料屬性設置為銅，球形空氣域的材料屬性為空氣；在腦出血MIT中常用的激勵頻率在1～10 MHz[16]，將激勵頻率設置為1 MHz，線圈中激勵電流為1 A，并按表1設置模型中各頭部組織的電磁特性參數(數據來源文獻[24])。

表1 1 MHz下的腦組織電磁特性

在腦模型中分別設置了3個不同位置的腦出血，分別位于空間位置A(0, 0.03, 0.075)處，位置B(-0.026, -0.04, 0.075)處及位置C(0.02, 0.02, 0.075)處，在每個位置分別設置出血量為24 ml, 14 ml及2 ml的腦出血。設置的腦出血分布情況如圖4所示。

4.2 仿真計算

4.2.1 相位差仿真計算

采用循環激勵，循環檢測的模式，即其中1個線圈為激勵，其余7個進行檢測，當遍歷8個線圈進行激勵后共得到56個檢測相位值。仿真計算圖4中位置A, B, C處3個出血量的相位差如圖5所示。

圖5 位置A,B,C處3個出血量的相位差

圖5(a)-圖5(c)分別表示位置A, B, C處3個出血量的相位差，每個圖中黑色、藍色、紅色虛線分別表示在該位置24 ml, 14 ml, 2 ml出血量情況下的仿真計算得到相位差測量值，由圖5(a)可以看出，檢測相位差數量級很小；腦出血體積的大小對線圈測量值影響很大，出血量越大檢測相位差的幅值變化越大，出血量越小檢測相位差的幅值變化越小，位置B, C處不同出血量的測量相位差規律與此類似。

4.2.2 靈敏度仿真計算

靈敏度矩陣是圖像重建的先驗信息[15]，常見的靈敏度矩陣求解方法有擾動法和互易定理法，采用基于互易定理的計算方法能大量減少計算時間[16]，靈敏度矩陣通過圖3中建立的腦出血MIT有限元模型仿真計算得到，假設有激勵-檢測線圈對a和b，分別對a和b施加激勵，另外一個線圈作檢測，此時可計算場域內各網格上電導率發生變化時的靈敏度，計算公式為

其中，NS是檢測線圈數量； Im(·) 代表虛部；Ωi是成像區域內網格單元；和分別是線圈a和線圈b施加激勵時檢測區域內的電場分布； (U0)a,b是空場時由激勵-檢測線圈ab計算得到的感應電壓。建立的模型是8線圈傳感器，因此檢測線圈數量NS=7。在參考頻率為1 MHz的情況下，通過將整個檢測區域設置為1 S/m的均勻電導率來計算靈敏度矩陣S。通過有限元仿真計算相關場量并根據式(10)計算靈敏度矩陣，計算6種典型線圈組合下的靈敏度如圖6所示。

圖6 不同激勵-檢測組合的靈敏度圖

圖6(a)-圖6(f)分別表示線圈1-線圈5、線圈4-線圈8、線圈2-線圈3、線圈5-線圈6、線圈5-線圈8、線圈3-線圈6組合時的靈敏度圖，每兩組分別代表相對、相鄰、相隔線圈激勵檢測情況下的典型靈敏度分布情況，其中不同的顏色反映該區域對電導率變化的敏感度不同，顏色對應數值的絕對值越大表示該區域靈敏度越高，顏色對應數值的絕對值越小表示該區域靈敏度越低；從圖6可以看出，在單線圈激勵單線圈檢測的情況下，越靠近激勵或檢測線圈的區域，靈敏度的絕對值越大，遠離該區域靈敏度絕對值逐漸減小，如圖6(a)中，線圈1-線圈5組合分別作激勵或檢測線圈時，越靠近線圈1和線圈5區域的靈敏度絕對值越大，該區域發生腦出血引起顱腦內電導率分布發生變化時，檢測線圈測量到的信號變化也越大，因此腦出血在靠近線圈的區域時成像效果更好，精度也更高。

5 實驗及結果分析

為了驗證所提改進NR算法的有效性，建立了真實顱腦和線圈的有限元模型，并設置了圖3中9種情況的腦出血，通過模型仿真計算得到9種腦出血情況下的線圈檢測相位差和靈敏度矩陣后，應用所改進的NR算法進行圖像重建實驗，并與NR算法[26]、Tikhonov算法[8]、CGLS算法[27]、Landweber算法[28]和優化迭代策略的NR算法[19]進行比較，采用改進的NR算法和5種對比算法分別對z=0.075平面進行斷層成像，圖像重建結果如圖7所示。實驗在處理器Intel i7-10875H，內存16 GB，顯卡Nvidia Ge-Force RTX 2 060, 6 GB顯存的電腦上進行。

在圖7(a)-圖7(c)、圖7(d)-圖7(f)、圖7(g)-圖7(i)分別對應位置A處、B處及C處出血量分別為24 ml, 14 ml, 2 ml共9種情況的腦出血圖像重建結果，其中的黑色虛線表示腦出血的實際輪廓。從圖7可以看出，6種成像算法均能重建出不同出血情況的腦出血位置，在同一位置下，隨著出血量的減小，與之對應區域的重建電導率幅值也減小，這是因為出血量減小，檢測的相位差信號也相對減小；Tikhonov, Landweber, CGLS, NR、優化迭代策略的NR等5種算法的重建圖像，出血的范圍都大于實際范圍，而使用所提改進NR算法的重建圖像，出血范圍與實際范圍更吻合，重建圖像中腦出血周圍的偽影比Tikhonov, Landweber, CGLS, NR、優化迭代策略的NR等5種算法少，且腦出血與背景的邊界更清晰，對出血量的判斷也更準確，可有效實現對2 ml腦出血的圖像重建。

為了評價不同算法的圖像重建性能，采用圖像相關系數(Correlation Coefficient, CC)、圖像重建誤差(Image Error, IE)、歸一化均方距離(Normalization Mean Square Distance criterion, NMSD)、圖像重建時間和迭代次數作為評價指標。其中圖像重建時間是算法開始到重建圖像所花費的時間，迭代次數是目標函數值小于設定閾值ε時已進行迭代計算的次數，各評價指標的計算公式為

其中，Δσtrue是真實電導率分布，Δσrec是重建的電導率分布，Δσtrue和 Δσrec分別是 Δσtrue和 Δσrec的平均值。利用CC和NMSD評價重建圖像的分辨率[11]，CC越大，NMSD越小，表示重建圖像分辨率更優；CC越大，IE和NMSD越小，表示圖像重建質量越高，分別計算9種出血情況下6種算法重建圖像的相關系數、重建誤差、歸一化均方距離、圖像重建時間和迭代次數，并分析了所提改進NR算法中投影算子P對收斂速度提高和成像質量改善的影響，得到的結果如表2、表3和表4所示。

表2 相關系數

表3 圖像誤差/歸一化均方距離

表4 圖像重建時間(s)/迭代次數

表2、表3、表4分別是重建圖像的相關系數、圖像重建誤差和歸一化均方距離、圖像重建時間和迭代次數，其中表的第1列是9種腦出血分布情況，第2列-第8列是9種出血情況下6種算法重建圖像以及無投影算子的改進NR算法的相關系數、重建誤差和歸一化均方距離、圖像重建時間和迭代次數；對比3個表中有無投影算子的改進NR算法，可以看出施加投影算子P后，算法的相關系數提高，重建誤差和歸一化均方距離均減小，且圖像重建時間和迭代次數均減少，平均為NR算法的1/2，表明通過投影算子P施加物理意義上的約束后提高了重建圖像的成像質量和收斂速度；由表2和表3可以看出，對于9種腦出血情況，采用所提改進NR算法重建圖像的相關系數均大于其他5種算法重建圖像的相關系數，歸一化均方距離均小于其他算法的歸一化均方距離，重建圖像的分辨率更優，同時其重建誤差也均小于其他算法的圖像重建誤差，表明所提改進NR算法重建圖像的質量更高；在同一出血位置，隨著出血量的減小，各算法重建圖像的相關系數減小，重建誤差增大，這是因為當出血量較小時，檢測到的相位信號變化也較小，圖像重建難度增大，同時對于2 ml微小出血的重建誤差較大，這是由于出血量較小檢測信號微弱以及檢測信號由各個組織及其相互間耦合的信號構成，可通過研究更有效的圖像重建算法、在圖像重建過程中引入腦組織的結構先驗信息以及腦出血MIT信號的有效分解來進一步降低重建誤差。結合圖7所給出的圖像重建結果，同一位置處出血量減少時，所提改進NR算法得到的腦出血范圍更符合實際腦出血范圍，其余5種算法得到的腦出血范圍雖然有所減小，但依然大于實際范圍，同時所提改進NR算法的相關系數更高、重建誤差和歸一化均方距離更低，并且能實現2 ml腦出血的圖像重建，表明所提改進NR算法求解得到的顱腦截面電導率分布及其數值解更符合實際，求解過程更加穩定。

由表4可以看出，對于9種腦出血情況，改進NR 算法重建圖像的時間比Tikhonov, Landweber,CGLS等3種算法較多，這是因為Tikhonov僅需一步成像，Landweber和CGLS在迭代過程中只需計算目標函數的1階導數[29,30]，而NR算法每步迭代都需要計算其2階導數，導致圖像重建時間較多。但所提改進NR算法相比NR和優化迭代策略的NR算法所需的重建時間大幅減少，成像時間平均只需NR算法的1/3。結合表2和表3分析，與Landweber,CGLS, NR和優化迭代策略的NR算法4種迭代圖像重建算法相比，所提改進NR算法迭代次數最少，重建圖像質量最高，表明所提改進NR算法每一步迭代的效率更高，收斂性更好且能夠得到更高質量的重建圖像。

6 結束語

本文提出一種用于腦出血MIT檢測的改進NR圖像重建算法，首先建立了一個包含頭皮、顱骨、腦脊液和腦實質4種組織的真實3維顱腦仿真模型，并仿真計算了9種腦出血情況下的線圈檢測值和靈敏度矩陣，采用改進的NR算法和其他5種圖像重建算法開展了圖像重建實驗。實驗結果表明，對于不同腦出血分布情況，所提改進NR算法的重建圖像相比于其他5種圖像重建算法在相關系數上提高了0.013～0.321，在重建誤差上減小了0.019～1.165，在歸一化均方距離上減小了0.013～1.060，成像時間平均只需NR算法的1/3，改進NR算法能實現較少的迭代次數和時間重建出更高質量的圖像，且能實現2 ml微小腦出血的圖像重建，表明所提改進NR算法的迭代效率更高、穩定性更好、重建圖像質量更高，由此可見，所提改進NR算法用于腦出血MIT是可行且有效的，有力推動了MIT技術的發展。

在本研究中，考慮到仿真模擬所需的計算資源，簡化了腦組織僅保留了主要的頭部組織，而真實腦組織更為復雜，同時改進NR算法重建圖像中腦出血和背景之間仍存在一定的偽影，下一步研究將構建更為精細的腦組織模型，并搭建腦出血MIT實驗平臺，利用瓊脂塊或3D打印技術制作不同腦出血情況下的真實3維顱腦模型，采集相關實驗數據，進一步驗證所提改進NR算法有效性，并開展應用研究。