基于改進對比源反演算法的腦卒中微波檢測成像技術(shù)研究

劉珩，王若璇

（北京理工大學(xué) 信息與電子學(xué)院, 北京 100081）

腦卒中是一種具有高病發(fā)率、高致殘率、高死亡率的急性腦血管疾病，是導(dǎo)致中國人口死亡的主要疾病之一[1]，據(jù)調(diào)研2018 年中國腦血管疾病已造成157 萬人死亡[2]. 目前國內(nèi)對腦卒中檢測的方式主要有CT 掃描、X 射線攝影、MRI 成像等[3],然而這些檢測手段存在著檢測設(shè)備龐大且昂貴、檢測時間長、多存在輻射危害等缺陷，不僅為患者帶來了巨大的經(jīng)濟負擔，甚至可能因為錯過救治的“黃金3 小時”從而給患者帶來不可逆轉(zhuǎn)的傷害[4]. 因此，快速、準確、有效的腦卒中檢測方式和血塊定位方法是救治腦卒中的關(guān)鍵.

針對腦卒中檢測現(xiàn)狀并結(jié)合微波技術(shù)特點，國內(nèi)外學(xué)者開始研究利用微波檢測腦卒中的新技術(shù).微波檢測技術(shù)依賴于不同組織之間介電對比度的存在[5]，具有使用安全，成像用時短，設(shè)備造價低廉且便于攜帶等優(yōu)點，可應(yīng)用于乳腺癌、腦卒中等病癥的檢測中. 近年來國內(nèi)外也出現(xiàn)了大量關(guān)于將微波成像技術(shù)應(yīng)用于腦卒中檢測的研究，所涉及到的微波成像技術(shù)包括Born 近似迭代法[6]、微波共焦成像法[7]、對比源反演成像法[8-9]、模式識別算法[3，10-11]等，關(guān)于微波檢測成像系統(tǒng)搭建的研究也層出不窮[12-14].2016 年，ISMAIL 等[9]首次提出將基于對比源反演(contrast source inversion, CSI)方法的微波成像技術(shù)應(yīng)用于對大腦中血塊的成像檢測中，并取得了良好的成像效果.

然而由于CSI 算法收斂速度較慢，導(dǎo)致成像的時間較長. 為了提高成像的速度，本文提出一種改進的對比源反演算法-FCG-CSI 算法，用于出血性腦卒中的微波檢測成像，利用快速共軛梯度法(fast conjugate gradient, FCG)替代常規(guī)共軛梯度法進行目標函數(shù)的優(yōu)化，從而加快CSI 算法的收斂. 建立腦卒中微波檢測成像系統(tǒng)，仿真實驗包括利用多物理場仿真軟件COMSOL 進行腦部模型的建立、設(shè)計天線進行微波信號的收發(fā)的模擬并采集相應(yīng)的微波數(shù)據(jù)、利用FCG-CSI 算法對收集到的數(shù)據(jù)進行反演成像.成像結(jié)果表明，F(xiàn)CG-CSI 算法相比于CSI 算法成像時間更短、成像精確度也更高.

1 腦卒中微波檢測成像整體框架

腦卒中微波檢測系統(tǒng)采用電磁逆散射的成像方法，如圖1 所示，系統(tǒng)可分為電磁散射正問題與電磁散射逆問題兩部分. 其中正問題是在已知入射場和對比度函數(shù)的情況下，計算求解空間中的總電場與散射場分布[15]. COMSOL 軟件是一款基于有限元法進行多物理場仿真計算的專業(yè)軟件，本文將利用COMSOL 仿真軟件進行電磁散射正問題的仿真求解，計算得到腦部周圍的電場分布，并對微波數(shù)據(jù)進行采集、處理與存儲；電磁散射逆問題根據(jù)入射場與計算得到的散射場，并通過重建算法反演人腦內(nèi)部結(jié)構(gòu). 由于逆散射問題是非線性且不適定的[15]，求解難度很大，因此重建算法的選擇對成像結(jié)果至關(guān)重要.本文利用一種改進的對比源反演算法-FCG-CSI算法對微波數(shù)據(jù)進行反演成像.

圖1 腦卒中微波檢測成像系統(tǒng)結(jié)構(gòu)Fig. 1 Structure of microwave imaging system for stroke detection

2 成像算法

2.1 CSI 算法基本理論

考慮如圖2 所示的電磁逆散射模型. 假設(shè)在一個復(fù)介電常數(shù)已知的均勻背景介質(zhì)中，有一形狀、大小、復(fù)介電常數(shù)均未知的散射體存在于有界且單連通的目標區(qū)域D. 發(fā)射天線與接收天線按如圖2 所示的圓形軌道分布在目標區(qū)域D周圍，發(fā)射天線數(shù)量為j，構(gòu)成測量區(qū)域S. 本文考慮在TM 極化的條件下，此時散射場可表示為[16]

圖2 電磁逆散射問題示意圖Fig. 2 Schematic diagram of electromagnetic inverse scattering problem

在以上積分方程中，式(1)被稱為“數(shù)據(jù)”方程，式(3)被稱為“目標”方程[17]. 為了可以在算法中更加簡潔且緊湊地表示“數(shù)據(jù)”與“目標”方程，引入對比源的概念，其為對比度與總電場的乘積，定義為w(r′)=χ(r′)E(r′)，這樣便可將“數(shù)據(jù)”方程和“目標”方程簡寫為積分算子的形式，即

對比源反演(CSI)算法采用常規(guī)共軛梯度法（即P-R 共軛梯度法）對代價函數(shù)進行迭代優(yōu)化，從而不斷更新對比源與對比度.

乘法正則化對比源反演（multiplicative regularized contrast source inversion，MR-CSI）算法是在CSI 算法的基礎(chǔ)上加入乘法正則化項FDTV(χ)，其代價函數(shù)為[16]

2.2 FCG-CSI 成像算法

在腦卒中微波檢測成像技術(shù)中，MR-CSI 算法具有良好的成像效果，但由于常規(guī)CSI 算法應(yīng)用的PR共軛梯度法收斂速度較慢，從而導(dǎo)致反演成像效率較低，耗時較長. 為了加快MR-CSI 算法的收斂速度，進一步提高算法的反演成像效率，本文在MR-CSI 算法的基礎(chǔ)上，利用快速共軛梯度法（FCG）來替代常規(guī)共軛梯度法來進行代價函數(shù)的優(yōu)化，由此來加速代價函數(shù)的收斂.

在常規(guī)MR-CSI 算法中，采用的是常規(guī)共軛梯度法，即PR（Polak-Ribiere）共軛梯度法來更新對比源.而在FCG-CSI 算法中，首先需將對比源wj,n-1進行一定的處理，以下用參數(shù)qj,n-1來表示經(jīng)處理后的對比源wj,n-1. FCG-CSI 算法中對比源的迭代公式為[18]

為了避免當代價函數(shù)曲線趨于收斂時出現(xiàn)上下抖動的現(xiàn)象，于是選擇在迭代趨于收斂時改為利用常規(guī)共軛梯度法進行迭代優(yōu)化.

綜上所述，當發(fā)射天線數(shù)量為j，頻率為k時，F(xiàn)CG-CSI 算法的步驟如下.

Step1: 輸入散射場數(shù)據(jù)fj,k與入射場數(shù)據(jù).

Step2: 初始化對比源wj,k,0與對比度χk,0.

ηS和ηD為權(quán)重因子，S表示測量區(qū)域，D表示目標區(qū)域，其表達式分別為

Step6: 當代價函數(shù)達到設(shè)定的誤差值或迭代次數(shù)達到預(yù)設(shè)的最大迭代次數(shù)時，輸出對比度矩陣χn并進行成像，算法結(jié)束，否則返回Step3.

3 仿真系統(tǒng)構(gòu)建與成像實現(xiàn)

3.1 腦部模型建立與微波數(shù)據(jù)獲取

運用多物理場仿真軟件COMSOL 仿真人體腦部模型與發(fā)射天線，組成腦卒中微波檢測成像系統(tǒng)的電磁散射正問題部分，進行微波數(shù)據(jù)的采集. 如圖3所示，人腦模型采用均勻簡易模型，即將腦部組織的相對介電常數(shù)統(tǒng)一為腦灰質(zhì)和腦白質(zhì)兩種組織的相對介電常數(shù)的平均值，腦部的半徑設(shè)為9 cm，內(nèi)部血塊的半徑設(shè)為1 cm.

圖3 微波腦部成像系統(tǒng)仿真模型Fig. 3 Simulation model of microwave brain imaging system

腦部各組織的電磁特性可通過DebyeⅡ 模型分析計算得到. DebyeⅡ 公式可表示為[19]

當入射微波信號的頻率在1～4 GHz 變化時，人體腦部組織各項電磁模型參數(shù)如表1 所示，將各項參數(shù)代入DebyeⅡ模型公式中，可得到1～4 GHz 頻率范圍內(nèi)腦組織與血塊的平均相對介電常數(shù)與平均電導(dǎo)率. 最終經(jīng)計算選取人腦均勻簡易模型的相對介電常數(shù)為40，電導(dǎo)率選取為0.9 S/m；血塊的相對介電常數(shù)選取為58，電導(dǎo)率選取為2.3 S/m.

表1 腦部各組織Debye Ⅱ模型參數(shù)[19-21]Tab. 1 Debye Ⅱ model parameters of brain tissues

仿真系統(tǒng)模型的最外層設(shè)置完美匹配層（PML層），以確保散射體產(chǎn)生的散射場幾乎完全被吸收，而不會有剩余部分反射到模型的外部邊界上. PML層與人腦模型之間為匹配介質(zhì)，匹配介質(zhì)的介電常數(shù)應(yīng)接近于腦部皮膚的介電常數(shù)，以減少天線和腦部之間的信號反射，從而確保電磁能量能夠進入大腦內(nèi)部，在這里選取匹配介質(zhì)的相對介電常數(shù)為42，電導(dǎo)率選取為0.9 S/m.

發(fā)射天線采用微帶貼片天線，其由一層薄金屬、一個矩形FR4 介質(zhì)塊和一個接地面組成. 微帶饋線、天線輻射器和接地面被建模為理想電導(dǎo)體表面，由50 Ω 的集總端口饋電來代表電源饋電. 由于從邊緣向貼片天線饋電會導(dǎo)致輸入阻抗非常高，進而導(dǎo)致不合需要的阻抗失配，因此本文中饋電方式選用嵌入式饋電，在貼片邊緣與中心之間尋找一最佳饋電點，以改善50 Ω 饋電點與天線之間的匹配情況. 切口區(qū)域的寬度W選為盡可能大，既使得天線與微帶之間的耦合最小又不會顯著影響天線的特性，所選取的饋電線長度L使反射功率S11最小，其最佳尺寸可通過參數(shù)化掃描與優(yōu)化得到. 經(jīng)設(shè)計與優(yōu)化最終得到天線結(jié)構(gòu)與其S參數(shù)曲線如圖4 所示，該天線的工作中心頻率為2 GHz，可通過設(shè)計與更改貼片的尺寸大小來設(shè)計滿足工作頻率與天線增益要求的天線.

圖4 天線示意圖Fig. 4 Diagram of antenna

在距離腦部模型中心10 cm 的圓形軌道上均勻放置16 個收發(fā)天線，每兩個天線的夾角為22.5°，如圖5 所示收發(fā)天線按照逆時針的方向從1～16 依次編號. 在數(shù)據(jù)收集時，所有天線依次單獨發(fā)射TM 波，COMSOL 軟件將根據(jù)數(shù)據(jù)方程和狀態(tài)方程，利用有限元法計算電磁場的時空分布，然后所有16 根天線進行微波數(shù)據(jù)的采集接收. 實驗需要測得區(qū)域內(nèi)存在散射體和不存在散射體兩種情況下的電場值，存在散射體時測得的微波數(shù)據(jù)為總電場值，散射體不存在時得到的微波數(shù)據(jù)為背景電場值，最終獲得21 616個微波數(shù)據(jù)，將這些數(shù)據(jù)導(dǎo)入MATLAB 中進行處理與反演成像.

圖5 收發(fā)天線分布Fig. 5 Distribution of transmit and receive antennas

3.2 成像結(jié)果與分析

腦卒中微波檢測系統(tǒng)的數(shù)據(jù)處理端使用MATLAB 進行運算得到反演成像結(jié)果，將獲得的微波數(shù)據(jù)導(dǎo)入MATLAB 中后，首先利用基于MATLAB 的FDFD 軟件包“MaxwellFDFD”，根據(jù)每一個發(fā)射天線對應(yīng)的16 個接收數(shù)據(jù)，計算得到該發(fā)射天線對應(yīng)的整個區(qū)域的散射電場分布矩陣與背景電場分布矩陣，其中散射電場可根據(jù)公式Etot=Einc+Esct，由總電場與背景電場相減得到. 然后將散射場矩陣與背景電場（入射場）矩陣作為FCG-CSI 算法的輸入進行反演成像.

為了量化成像算法的性能，利用 Δr指標來分析成像結(jié)果的位置誤差，令Lreal為真實散射體的幾何中心，Lmax為成像區(qū)域能量單元所在位置，S為目標散射體區(qū)域，Z為除去散射體區(qū)域外的成像區(qū)域，則用Δr來表示成像目標位置的精準度， Δr值越小，定位越準確. 位置誤差 Δr定義為

圖6 成像結(jié)果Fig. 6 Imaging results

經(jīng)計算得到分別在頻率為2 GHz、2.38 GHz 時，通過CSI 算法與FCG-CSI 算法得到的成像結(jié)果中，預(yù)測血塊的中心位置與實際模型中血塊中心位置的偏差如表2 所示. 可以看出在兩種頻率下，F(xiàn)CG-CSI算法相對于常規(guī)CSI 算法的 Δr值均更小，成像結(jié)果更加精確.

表2 成像結(jié)果位置誤差Tab. 2 Position error of imaging results

圖7 顯示了隨著迭代次數(shù)的增加，利用FCGCSI 算法和CSI 算法進行反演成像時目標函數(shù)的收斂曲線，可以看出在迭代初期，F(xiàn)CG-CSI 算法比CSI算法收斂速度更快. 目標函數(shù)反映的是通過迭代計算得到的“數(shù)據(jù)”與“目標”方程的值（即計算得到的散射場與總場的值）與實際測量的散射場與總場的值之間的誤差，經(jīng)過多次實驗發(fā)現(xiàn)，當此誤差值小于0.2 時基本可以完成準確成像，再往下迭代對成像結(jié)果的影響不大. 目標函數(shù)達到0.2 以下常規(guī)CSI 算法需要7 次迭代，而FCG-CSI 算法僅需要5 次迭代，因此FCG-CSI 算法的成像速度相比CSI 算法更快.

圖7 迭代次數(shù)-代價函數(shù)曲線圖Fig. 7 Curve of iterations number-cost function

為了進一步提高成像的精度，可以利用多個頻率下得到的微波數(shù)據(jù)進行反演成像，其方法是在FCG-CSI 算法步驟的式(11)(13)(15)(17)(21)中·運算的基礎(chǔ)上添加運算符，進行不同頻率微波數(shù)據(jù)的疊加. 由于結(jié)合了多頻微波數(shù)據(jù)信息，因此理論上可得到更精準的成像結(jié)果. 圖8 顯示了FCG-CSI算法在1.3 GHz、1.5 GHz、2 GHz、2.38 GHz 頻率下，結(jié)合多頻微波數(shù)據(jù)的反演成像結(jié)果. 此時的位置偏差值為Δr=5.385 2 mm. 可以看出此時的成像位置偏差均小于各單獨頻率下的成像結(jié)果，這說明利用多頻率數(shù)據(jù)可以得到比單頻率更精準的成像效果.

圖8 多頻數(shù)據(jù)FCG-CSI 算法成像結(jié)果Fig. 8 Imaging results of FCG-CSI algorithm for multi-frequency data

4 結(jié)束語

本文將對比源反演算法應(yīng)用于出血性腦卒中的微波無損檢測成像中，并在對比源反演算法的基礎(chǔ)上加以改進，利用快速共軛梯度法來替代常規(guī)共軛梯度法來進行代價函數(shù)的優(yōu)化，從而加速代價函數(shù)的收斂，提高反演成像的效率. 使用COMSOL 仿真軟件進行腦部的建模與微波信號收發(fā)過程的仿真，然后將采集的微波數(shù)據(jù)導(dǎo)入MATLAB 軟件中進行反演成像. 從成像結(jié)果來看，反演成像所使用的FCGCSI 算法相對于常規(guī)的CSI 算法，不僅收斂速度更快，成像結(jié)果也更加精確. 在利用多頻數(shù)據(jù)進行FCGCSI 算法反演成像時，成像的位置誤差可達到0.54 cm.

本文的研究均是在較為理想的仿真環(huán)境中完成，并未考慮實際測量環(huán)境中噪聲的影響. 在后續(xù)研究中，將考慮到噪聲對微波檢測成像的影響并設(shè)法改善成像算法的抗噪聲性能. 與此同時，在接下來的研究中也將進一步細化腦部模型，使仿真更加貼近于真實檢測.