劉 俐 李 倩 何 為* 徐 征

1(重慶市腫瘤研究所，重慶 400030)2(重慶大學(xué)電氣工程學(xué)院 輸配電裝備及系統(tǒng)安全與新技術(shù)國家重點(diǎn)實(shí)驗(yàn)室，重慶 400044)

一種均勻激勵磁場磁感應(yīng)成像的改進(jìn)反投影算法

劉 俐1李 倩2何 為2*徐 征2

1(重慶市腫瘤研究所，重慶 400030)2(重慶大學(xué)電氣工程學(xué)院 輸配電裝備及系統(tǒng)安全與新技術(shù)國家重點(diǎn)實(shí)驗(yàn)室，重慶 400044)

介紹一種均勻激勵磁場磁感應(yīng)成像的改進(jìn)反投影算法。根據(jù)畢奧-薩伐爾定律和靈敏度的分布，提出以靈敏度為零的點(diǎn)作為分界點(diǎn)，采用不同方式的加權(quán)系數(shù)，對反投影算法進(jìn)行修正。靈敏度分布的仿真實(shí)驗(yàn)表明，當(dāng)檢測線圈與成像區(qū)間的剖分單元之間的距離r等于0.06 m時(shí)，檢測線圈上的測量電壓的變化對距離r的變化的靈敏度為0，此時(shí)剖分單元的位置標(biāo)記為P0。當(dāng)檢測線圈與剖分單元的距離r小于0.06 m時(shí)，反投影算法的加權(quán)系數(shù)取為r2；反之r大于等于0.06 m時(shí)，加權(quán)系數(shù)均為P0到檢測線圈的距離r0的平方。不同位置的目標(biāo)的實(shí)驗(yàn)結(jié)果表明，當(dāng)異物位于成像區(qū)間的中心和偏離中心0.015 m時(shí)，都能準(zhǔn)確反映異物的位置，但是當(dāng)異物偏離中心0.03 m時(shí)，得到的成像結(jié)果比實(shí)際模型更偏離中心。不同大小的目標(biāo)的實(shí)驗(yàn)結(jié)果表明，與實(shí)際模型中異物尺寸相比，直徑較小的異物的成像基本能夠反映出模型中實(shí)際異物的大小，而直徑較大的異物的成像在尺寸上和位置上均難以體現(xiàn)模型中實(shí)際異物的真實(shí)形態(tài)。本研究提出的采用不同方式的加權(quán)系數(shù)的改進(jìn)反投影算法，能一定程度上消除傳統(tǒng)反投影算法中目標(biāo)成像拖影范圍較大的痼疾，對異物的定位較為準(zhǔn)確，是一種在磁感應(yīng)成像研究中值得推廣的新方法。

磁感應(yīng)成像；反投影算法；亥姆霍茲線圈

引言

磁感應(yīng)成像[1](magnetic induction tomography，MIT)是一種無創(chuàng)的非接觸功能型成像技術(shù)，采用時(shí)變磁場作為激勵源，以避免接觸阻抗的對成像造成的誤差。由于磁場具有穿透高電阻率物質(zhì)的特性[2]，能夠透過高電阻率的顱骨直達(dá)病變區(qū)域；磁場還具有高度的聚焦性，對成像區(qū)域中心位置有較強(qiáng)敏感性。因此，磁感應(yīng)成像技術(shù)更適合應(yīng)用于頭部疾病的成像，是目前生物阻抗成像技術(shù)領(lǐng)域的熱門研究課題。

磁感應(yīng)成像逆問題是指在已知激勵源和邊界感應(yīng)電壓的情況下，重構(gòu)出模型內(nèi)電導(dǎo)率的分布。這是一個高病態(tài)的非線性問題[3-4]，求解過程復(fù)雜，并且由于顱內(nèi)組織的低電導(dǎo)率和弱對比度，建立有效的重構(gòu)算法相對困難，國內(nèi)外的研究人員對其重構(gòu)算法進(jìn)行了許多研究。文獻(xiàn)[5]運(yùn)用修正的Newton-Raphson(N-R)圖像重建算法實(shí)現(xiàn)了單目標(biāo)和多目標(biāo)的圖像重建。文獻(xiàn)[6]在磁感應(yīng)成像的重構(gòu)算法中做了大量工作,比如Tikhonov正則化算法, 改進(jìn)的高斯牛頓一步動態(tài)重建算法。Scharfetter等將反投影算法運(yùn)用于磁感應(yīng)成像中[7]。運(yùn)用靈敏度矩陣算法的一個重要問題是求解靈敏度矩陣的廣義逆矩陣，Morris采用截?cái)嗥娈愔捣纸夥ㄇ蠼忪`敏度矩陣的廣義逆矩陣[8]。Casanova等的研究表明，采用Tikhonov正則化能夠提高截?cái)嗥娈愔捣纸馑惴ㄐ阅躘9]。呂等采用補(bǔ)償原理計(jì)算靈敏度矩陣[10]。以上的算法均涉及到大量的矩陣運(yùn)算，成像速度比較慢，而對于臨床進(jìn)行連續(xù)實(shí)時(shí)監(jiān)護(hù)顱內(nèi)組織病變，需要選用一種速度快的重構(gòu)算法。本研究采用改進(jìn)的濾波反投影算法，引入混合加權(quán)法。本研究針對傳統(tǒng)磁感應(yīng)成像系統(tǒng)激勵磁場的不均勻性對成像精度造成的影響，提出采用能夠產(chǎn)生均勻磁場的亥姆霍茲線圈作為激勵線圈。在此基礎(chǔ)上運(yùn)用反投影算法對單目標(biāo)模型進(jìn)行圖像重構(gòu)，并提出采用混合加權(quán)系數(shù)對反投影算法進(jìn)行修正，以消除成像偽跡。

1 均勻激勵磁場的MIT模型

目前已有的磁感應(yīng)成像系統(tǒng)通常采用單邊螺旋小線圈作為激勵線圈[11-12], 此種激勵方式產(chǎn)生的磁場在空間中的分布是發(fā)散和不均勻的，所以在計(jì)算逆問題時(shí)，需要對激勵磁場加復(fù)雜的權(quán)重，運(yùn)算更復(fù)雜，耗時(shí)多，影響成像速度。同時(shí)由于傳統(tǒng)的激勵磁場的深度比較淺，方向不一致，如圖1所示，這使相同容積和量值的電導(dǎo)率單元在成像區(qū)域中的位置不同，受到的激勵磁場大小不同，對圖像重構(gòu)造成影響。

為克服上述問題，采用能夠產(chǎn)生均勻磁場的亥姆霍茲線圈作為激勵線圈，其結(jié)構(gòu)如圖2所示，它由一對半徑和線圈匝數(shù)都相同的線圈組成，這兩個線圈彼此平行共軸放置，線圈間的中心距離等于它們的半徑。當(dāng)它們通入大小相等、方向相同的電流I時(shí)，這兩個載流線圈的總磁場在軸的中點(diǎn)附近較大范圍內(nèi)是均勻的，如圖3所示。

由此所提出的基于亥姆霍茲線圈激勵的均勻激勵磁場MIT系統(tǒng)，其結(jié)構(gòu)如圖4所示。其中e1和e2為一組由亥母赫茲線圈構(gòu)成的激勵線圈，d1～d8為相同的8個檢測線圈，被測物放置在均勻磁場區(qū)域。亥姆霍茲線圈在空間中產(chǎn)生均勻的激勵磁場B0，方向近似為直線，當(dāng)激勵磁場穿過導(dǎo)電物體時(shí)，將在物體內(nèi)部產(chǎn)生渦流，從而產(chǎn)生二次磁場ΔB擾動主磁場，在檢測線圈上將檢測到主磁場和二次磁場的合成場B0+ΔB產(chǎn)生的感應(yīng)電壓。由于不同的電導(dǎo)率分布產(chǎn)生的渦流以及ΔB的大小不同，因此檢測到的感應(yīng)電壓也不一樣。對檢測電壓數(shù)據(jù)進(jìn)行反投影，就能重構(gòu)出電導(dǎo)率的分布。

假設(shè)在檢測線圈和激勵線圈之間存在一條由磁力線構(gòu)成的投影帶，并且投影帶外的電導(dǎo)率變化對投影值沒有貢獻(xiàn)，再將檢測線圈上的電壓沿著投影帶的方向疊加回去。通過采用亥母赫茲線圈作為激勵線圈，使傳統(tǒng)激勵方式下彎曲的磁力線變成直線，由此投影路徑變得簡單，反投影算法更容易實(shí)現(xiàn)。

MIT基本測量原理如圖5所示。每次測量m個檢測線圈的電壓值Um(m為檢測線圈個數(shù))。再將檢測系統(tǒng)以亥姆霍茲線圈中點(diǎn)O為中心旋轉(zhuǎn)一定角度ΔΦ(如圖5中由角度1旋轉(zhuǎn)到角度2，即位置1旋轉(zhuǎn)到位置2)，得到另一組電壓數(shù)據(jù)。如此重復(fù)，直至旋轉(zhuǎn)n次，使nΔΦ=360°為止，即得到m×n維的電壓值矩陣。根據(jù)上述的測量方式，分別對存在異物的情況和電導(dǎo)率均勻分布時(shí)的模型進(jìn)行測量，分別得到電壓矩陣Up和Uu，這兩者之差即為異物擾動電壓Ug。將Ug沿激勵線圈與檢測線圈之間直線磁力線構(gòu)成的投影帶進(jìn)行反投影，就是令該區(qū)域的像素值為Ug。將每組測量擾動電壓沿其所對投影帶反投影，就可得到由磁感應(yīng)成像反投影算法得到的電導(dǎo)率重構(gòu)圖。反投影算法的計(jì)算過程表示為

(1)

式中，Cp[e]表示第e個剖分單元的電導(dǎo)率，E表示剖分單元的總數(shù)量。

2 改進(jìn)反投影算法-混合加權(quán)法消除圖像偽跡

根據(jù)畢奧-薩伐爾定律，源電流回路在空間某點(diǎn)產(chǎn)生的磁感應(yīng)強(qiáng)度與該點(diǎn)到源電流回路的距離的平方成反比。為了消除這種影響，可以選用有限元剖分單元到檢測線圈的距離r的平方作為加權(quán)系數(shù)，對式(1)提出的反投影公式加以修正，即

(2)

由于MIT成像技術(shù)中，電導(dǎo)率異常區(qū)域越靠近檢測線圈，其感應(yīng)電壓越大，為了研究異物位置對檢測線圈電壓的影響，分析了該模型的靈敏度分布。

靈敏度原理的主要思想是基于Geselowitz的阻抗電場理論，利用有限元方法將場域等效為一電阻網(wǎng)絡(luò)模型，根據(jù)這一網(wǎng)絡(luò)模型建立總體傳輸阻抗變化與單個單元中電阻變化之間的對應(yīng)關(guān)系，進(jìn)而利用這一對應(yīng)關(guān)系來修正阻抗分布[13,14]，靈敏度矩陣中靈敏度系數(shù)表示為

(3)

式中，Sij表示的是第i對電極對第j個單元的電導(dǎo)率的靈敏度系數(shù)，如圖6所示，區(qū)域Ω的電導(dǎo)率為σ，電極對(A,B)和(C,D)貼放在區(qū)域邊界?Ω上。設(shè)U表示激勵電流IAB加到電極對(A,B)上時(shí)區(qū)域的電壓分布，UCD表示此時(shí)電極對(C,D)測量到的電壓。V表示激勵電流ICD加到電極對(C,D) 上時(shí)區(qū)域的電壓分布，VAB表示此時(shí)電極對(A,B)測量到的電壓。根據(jù)式(3)，可仿真得到檢測線圈的靈敏度分布圖。

3 實(shí)驗(yàn)仿真

4 結(jié)果

4.1檢測線圈的靈敏度分布

根據(jù)式(3)，仿真得到檢測線圈的靈敏度分布圖如圖8所示，其靈敏度分布圖是渦流電流分布密度的相對值。圖8中(a)和(b)為不同檢測線圈的被測區(qū)域靈敏度的分布，可以看出在靠近線圈位置的靈敏度值最大，隨著電導(dǎo)率單元到檢測線圈距離的增加，其值迅速衰減。這表明在靠近線圈的位置被測物的電導(dǎo)率變化對檢測電壓影響大，遠(yuǎn)離線圈的電導(dǎo)率變化影響很小。為了分析靈敏度隨距離的變化情況，取每個檢測線圈軸線所對的直線上的靈敏度值，繪制其變化曲線，如圖9所示。從圖9中可以看出靈敏度值呈指數(shù)衰減，到距檢測線圈0.06 m(點(diǎn)P0處)外的靈敏度值為零。這表明P0外的電導(dǎo)率擾動對檢測線圈的電壓無影響。因此，選用混合加權(quán)的方法，即電導(dǎo)率單元到檢測線圈的距離小于點(diǎn)P0到檢測線圈的距離時(shí)，采用平方加權(quán)法；大于P0到檢測線圈的距離時(shí)的加權(quán)值均等于P0到檢測線圈的距離r0的平方，此處，即

(4)

4.2反投影算法仿真

測量角度下的擾動電壓數(shù)據(jù)按上述方法全部進(jìn)行反投影后就能得到仿真模型的重構(gòu)圖像，結(jié)果如圖10所示。其中靠近圓心的深色區(qū)域正好是模型中電導(dǎo)率異常區(qū)域所在的位置，表明該算法能夠準(zhǔn)確定位異物的位置。但是利用這種方法得到的重建圖像邊界很不清晰，存在星狀偽跡。這是因?yàn)榉赐队八惴ǖ膶?shí)質(zhì)是將n次投影數(shù)據(jù)均勻地疊加到每個剖分單元上去，使原本為零的點(diǎn)值變得不為零。

4.3改進(jìn)后反投影算法仿真結(jié)果

利用混合加權(quán)的反投影算法，并通過對不同模型的仿真結(jié)果來考察算法的性能。對不同位置和不同體積的異物，按式(4)進(jìn)行加權(quán)的圖像重構(gòu)，其成像結(jié)果如圖11和圖12所示，與圖10中反投影算法得到的結(jié)果對比可以看出，混合加權(quán)的反投影算法消除了部分偽跡，成像效果明顯改善。從圖11中可以看出，與原模型相比，經(jīng)過混合加權(quán)后的反投影算法的成像結(jié)果中異物的中心位置均與原模型設(shè)置的中心位置其坐標(biāo)位置一致。從圖12中可以看出，與原模型相比，混合加權(quán)算法能夠地對異物進(jìn)行定位，成像異物的大小順序也與實(shí)際模型的大小順序的基本一致。

用混合加權(quán)的反投影算法對3種不同位置的異物模型重構(gòu)所得到的圖像，與實(shí)際模型相比，成像結(jié)果中異物的大小均與實(shí)際模型設(shè)置的大小相符合。但是在異物的定位方面，當(dāng)異物位于中心和偏離中心0.015 m的情況下都能準(zhǔn)確反映異物的位置，但是當(dāng)異物偏離中心0.03 m的情況下，得到的結(jié)果比實(shí)際模型更偏離中心。

5 討論與結(jié)論

本研究提出了一種利用混合加權(quán)改進(jìn)磁感應(yīng)成像反投影算法，對比了改進(jìn)前后的兩種成像結(jié)果。在選取不同加權(quán)系數(shù)時(shí)，主要根據(jù)靈敏度的分布來確定。當(dāng)檢測線圈與剖分單元的距離r小于0.06 m時(shí)，反投影算法的加權(quán)系數(shù)取為r2；反之r大于等于0.06 m時(shí)，加權(quán)系數(shù)均為P0到檢測線圈的距離r0的平方。不同位置的目標(biāo)的實(shí)驗(yàn)結(jié)果表明，當(dāng)異物位于成像區(qū)間的中心和偏離中心0.015 m時(shí)，都能準(zhǔn)確反映異物的位置，但是當(dāng)異物偏離中心0.03 m時(shí)，得到的成像結(jié)果比實(shí)際模型更偏離中心。不同大小的目標(biāo)的實(shí)驗(yàn)結(jié)果表明，與實(shí)際模型中異物尺寸相比，直徑較小的異物的成像基本能夠反映出模型中實(shí)際異物的大小，而直徑較大的異物的成像在尺寸上和位置上均難以體現(xiàn)模型中實(shí)際異物的真實(shí)形態(tài)。

同時(shí),實(shí)驗(yàn)結(jié)果也表明，由于異物的存在產(chǎn)生渦流并攜帶異物電導(dǎo)率信息，該渦流產(chǎn)生的二次場對于檢查線圈的貢獻(xiàn)呈非線性，從而當(dāng)物體在區(qū)域中部時(shí)異物重構(gòu)的結(jié)果相對發(fā)散，而當(dāng)異物靠近邊緣時(shí)，檢測線圈獲取的信號相關(guān)性較好，重構(gòu)圖像相對集中，總體來說，成像目標(biāo)大小為實(shí)際物體的2至3倍，并與背景與物體的電導(dǎo)率差異相關(guān)。在后續(xù)的工作中，將根據(jù)這一關(guān)系建立校準(zhǔn)體系，對重構(gòu)圖像的等位線進(jìn)一步劃分，以期使重構(gòu)結(jié)果與物體尺寸更切合。此外，研究中假定投影帶外的電導(dǎo)率變化對感應(yīng)線圈沒有貢獻(xiàn)，這種簡化假定與實(shí)際并不相符，而采用上述受影響的感應(yīng)線圈結(jié)果進(jìn)行重構(gòu)分析必然導(dǎo)致范圍擴(kuò)大，要解決這該問題可以一方面增大線圈分布密度，另一方面將每個感應(yīng)線圈上的信號分解為其本身與周邊線圈結(jié)果的加權(quán)所得，對結(jié)果予以改善。

本研究提出的采用不同方式的加權(quán)系數(shù)的改進(jìn)反投影算法，能一定程度上消除傳統(tǒng)反投影算法中目標(biāo)成像拖影范圍較大的痼疾，對異物的定位較為準(zhǔn)確，是一種在磁感應(yīng)成像研究中值得推廣的新方法。

[1] Griffiths H, Stewart WR, Gough W. Magnetic induction tomography: a measuring system for biological tissues [J]. Ann NY Acad Sci, 1999,873:335-345.

[2] Netz J, Fornet E, Elaagecnannf S. Contactless impedance measurement by magnetic induction apossible method for investigation of brain impedance [J]. Physiol Meas, 1993,14:463-471

[3] Merwa R, Hollaus K, Brunner P,etal. Solution of the inverse problem of magnetic induction tomography (MIT) [J]. Physiol Meas, 2005,26: 241-250.

[4] Casanova R, da Silva AF, Borges AR,etal. MIT image reconstruction based on edge preserving regularization [J]. Physiol Meas, 2004,25:195-207.

[5] 王聰,劉銳崗,李燁,等. 一種用于磁感應(yīng)斷層成像的圖像重建算法 [J].儀器儀表學(xué)報(bào), 2008,29(10):2052-2057.

[6] 呂軼，王旭，金晶晶，等. 正則化一步動態(tài)重建算法在磁感應(yīng)成像中的應(yīng)用[J]. 電子學(xué)報(bào), 2011,39(12):2801-2806.

[7] Scharfetter H, Riu P, Populo M. Sensitivity maps for low-contrast perturbations within conducting background in magnetic induction tomography [J]. Physiol Meas, 2002,23(1): 195-202.

[8] Morris A, Griffiths H, Gough W. A numerical model for magnetic induction tomography measurements in biological tissues [J]. Physiol Meas, 2001,22:113-119.

[9] Casanova R, da Silva AF, Borges AR,etal. MIT image reconstruction based on edge preserving regularization [J]. Physiol Meas, 2004,25:195-207.

[10] 呂軼, 王旭, 楊丹,等. 一種用于磁感應(yīng)成像中靈敏度矩陣的計(jì)算方法[J].東北大學(xué)學(xué)報(bào)(自然科學(xué)版), 2011,32(5):618-621.

[11] Watson S, Williams RJ, Griffiths H,etal. Frequency downconversion and phase noise in MIT [J]. Physiol Meas, 2002,23:189-194.

[12] Scharfetter H, Lackner HK, Rosel J. Magnetic induction tomography: hardware for multi-frequency measurements in biological tissues [J]. Physiol Meas, 2001,22: 131-146.

[13] Geselowitz DB. An application of electrocardiographic lead theory to impedance plethysmography [J]. IEEE Transactions on Biomedical Engineering, 1971,18(1):38-40.

[14] 呂軼,王旭,金晶晶, 等.基于互易原理磁感應(yīng)成像中靈敏度矩陣的計(jì)算[J]. 儀器儀表學(xué)報(bào), 2012,33(3):616-623.

AnImprovedBack-ProjectionAlgorithmofMagneticInductionTomographywithUniformMagneticExcitation

LIU Li1LI Qian2HE Wei2*XU Zheng2

1(ChongqingCancerInstitute,Chongqing400030,China)2(StateKeyLaboratoryofPowerTransmissionEquipment&SystemSecurityandNewTechnology,TheElectricalEngineeringCollege,ChongqingUniversity,Chongqing400044,China)

This paper describes an improved back-projection algorithm of magnetic induction tomography with uniform magnetic excitation. According to Biot-Savart law and the distribution of the sensitivity, a new algorithm with mixed weight was proposed to correct back-projection algorithm. Considering the sensitivity of zero point as the dividing point, this algorithm used different weighting coefficients which were separated by the dividing point. According to the simulation experiment of sensitivity distribution, when the distancerfrom detection coil to the meshing element of imaging area was 0.06 m, the sensitivity value was zreo, and the position of this meshing element was markedP0. When the distancerfrom detection coil to meshing element was less than 0.06 m, the weighted coefficient wasr2; otherwise, the weighting coefficients were the square of the distancer0from the detection coil toP0point. The experimental results of different positions of foreign body showed that the back projection algorithm with mixed weight well reflected the size of the target object. When the foreign object was located in the center and 0.015 m off-center, the algorithm accurately reflected the position. However, when foreign object were 0.03 m off center, the location from the reconstructed image was more off-center than that of the actual model. The experimental results of different size of foreign body revealed that the reconstructed image of smaller diameter of the foreign object could accurately reflect the size of the foreign object, while that of larger diameter could not. The improved back projection algorithm can eliminate a wide range of smearing image errors caused by the conventional back-projection algorithm to a certain extent, and can also locate the foreign object more accurately. Therefore, the proposed algorithm is of significance in the study of magnetic induction image.

magnetic induction tomography; back projection algorithm; Helmholtz coil

10.3969/j.issn.0258-8021. 2014. 03.08

2013-03-30， 錄用日期:2014-05-06

國家自然科學(xué)基金(51107150)；中央高?；究蒲袠I(yè)務(wù)經(jīng)費(fèi)資助項(xiàng)目(CDJZR10150021)；“111”計(jì)劃項(xiàng)目(B08036)

R318.03

A

0258-8021(2014) 03-0313-07

*通信作者(Corresponding author)，E-mail: hewei@cqu.edu.cn