劉 俐 李 倩 何 為* 徐 征

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

一種均勻激勵磁場磁感應成像的改進反投影算法

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

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

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

磁感應成像；反投影算法；亥姆霍茲線圈

引言

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

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

1 均勻激勵磁場的MIT模型

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

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

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

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

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

(1)

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

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

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

(2)

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

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

(3)

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

3 實驗仿真

4 結(jié)果

4.1檢測線圈的靈敏度分布

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

(4)

4.2反投影算法仿真

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

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

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

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

5 討論與結(jié)論

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

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

本研究提出的采用不同方式的加權(quán)系數(shù)的改進反投影算法，能一定程度上消除傳統(tǒ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] 王聰,劉銳崗,李燁,等. 一種用于磁感應斷層成像的圖像重建算法 [J].儀器儀表學報, 2008,29(10):2052-2057.

[6] 呂軼，王旭，金晶晶，等. 正則化一步動態(tài)重建算法在磁感應成像中的應用[J]. 電子學報, 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] 呂軼, 王旭, 楊丹,等. 一種用于磁感應成像中靈敏度矩陣的計算方法[J].東北大學學報(自然科學版), 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] 呂軼,王旭,金晶晶, 等.基于互易原理磁感應成像中靈敏度矩陣的計算[J]. 儀器儀表學報, 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

國家自然科學基金(51107150)；中央高校基本科研業(yè)務經(jīng)費資助項目(CDJZR10150021)；“111”計劃項目(B08036)

R318.03

A

0258-8021(2014) 03-0313-07

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