馬虹霞 劉志朋 張順起 殷 濤

(中國醫(yī)學(xué)科學(xué)院 北京協(xié)和醫(yī)學(xué)院 生物醫(yī)學(xué)工程研究所，天津 300192)

引言

生物電阻抗成像技術(shù)［1－3］是近年發(fā)展起來的新型功能成像技術(shù)，它根據(jù)人體組織在不同生理和病理狀態(tài)下具有不同的電阻(導(dǎo))率，通過給人體施加安全驅(qū)動電流(電壓)，在體外測量響應(yīng)電壓(電流)信號，以重建人體內(nèi)部的電阻(導(dǎo))率分布和變化的圖像。生物電阻抗成像對組織電阻(導(dǎo))率的變化具有較高的檢測靈敏度和對比度，對疾病的早期診斷具有較明顯的優(yōu)勢，但空間分辨率較低的弱點限制了其在臨床上的深入應(yīng)用。為提高空間分辨率改善重建圖像的質(zhì)量，He等于2005年提出了一種新的組織功能成像技術(shù)—感應(yīng)式磁聲成像［4－5］(MAT-MI)。感應(yīng)式磁聲成像兼具電阻抗成像高對比度和超聲成像高空間分辨率的優(yōu)點，其基本原理是將目標物體置于靜態(tài)磁場中，施加磁脈沖激勵在目標物體中感應(yīng)出渦電流，渦電流與靜磁場作用產(chǎn)生洛倫茲力，帶電粒子受到洛倫茲力產(chǎn)生聲振動，由聲探頭檢測聲信號，通過圖像重建得到由聲信號攜帶的組織電導(dǎo)率分布與變化信息。

近年來感應(yīng)式磁聲成像作為一種新型功能成像技術(shù)備受關(guān)注，目前已有該方面的相關(guān)報道［6－15］，但是大多數(shù)研究都是建立電導(dǎo)率突變的規(guī)則模型并且假設(shè)激勵磁場空間分布均勻。由于生物組織結(jié)構(gòu)復(fù)雜，電磁場正問題中很難得到感應(yīng)渦流的解析解，并且組織電導(dǎo)率分布并非嚴格躍變，而是存在一定的過渡帶。因此，本研究從不規(guī)則的幾何模型出發(fā)，建立了具有一定尺寸的圓形單線圈以及正方體、偏心球、橢球和電導(dǎo)率連續(xù)分布的球模型，應(yīng)用有限元分析軟件 Comsol 3.5a，計算了具有不同電導(dǎo)率模型中的感應(yīng)渦流，利用 Matlab軟件工具根據(jù)有限元分析的數(shù)據(jù)，仿真計算了模型的聲源分布及逆問題重建所需的空間聲壓序列，應(yīng)用時間反演方法重建了能夠反映不同電導(dǎo)率邊界的聲源圖像信息。實際應(yīng)用中，由于客觀條件的限制，聲探頭掃描范圍有限，很難得到完備的測量數(shù)據(jù)，因此，本文最后還定性分析了探測器檢測范圍的有限性對重建聲源圖像的影響。

1 理論方程

感應(yīng)式磁聲成像是一種多物理場耦合的新型功能成像方法，涉及到電磁學(xué)、力學(xué)和聲學(xué)知識，研究MAT-MI可以分為正、逆問題兩個方面，正問題已知樣本電導(dǎo)率分布求解聲源及聲信號，逆問題是已知聲信號通過相應(yīng)的重建算法求解聲源進而得到樣本的電導(dǎo)率。

1.1 正問題

感應(yīng)式磁聲成像產(chǎn)生的聲信號滿足如下波動方程［4］

式中，cs是聲在組織中傳播的速度，p(r，t)是產(chǎn)生的聲波，J(r，t)是感應(yīng)渦流密度，r為被檢組織內(nèi)一點，Δ·(J(r，t)×B0)是聲振源。對于波動方程右側(cè)的聲振源可通過電磁場理論獲得，經(jīng)典麥克斯韋方程組為

本構(gòu)關(guān)系為

式(2c)表明磁通密度B是無散的，可表示成矢量磁位A的旋度，即有

將式(6)代入(2a)，有

亦即

在準靜態(tài)情況下，位移電流相對傳導(dǎo)電流可以忽略，將式(3)、(5)、(6)和式(9)代入式(2b)得到

由于樣本內(nèi)無外部電流源，并且生物電流遠小于感應(yīng)渦電流，根據(jù)電流連續(xù)性定理，有 Δ·J=0，即

因此，式(10)和式(11)構(gòu)成求解目標物體中感應(yīng)渦流的控制方程。當(dāng)激勵磁場脈寬足夠窄時，通過格林函數(shù)法求波動方程(1)得到聲壓的表達式為［4］

式中，rd為傳聲器檢測面上一點，Ω是聲傳感器檢測面。通過有限元法求解感應(yīng)渦流后，由式(12)便可獲得MAT-MI產(chǎn)生的聲信號，再利用此聲壓信號進行逆問題聲源圖像的重構(gòu)。

1.2 逆問題

假設(shè)激勵脈寬足夠窄時，感應(yīng)渦流可近似為J(r，t)=J(r)δ(t)，在時間域(－ ∞ ，0+)上對波動方程兩邊進行積分得到= Δ·(J×B0)，由聲場中的時間反演法可以得到聲源的重建公式為［4］

式中，Σ是傳聲器檢測面，r是被檢測組織內(nèi)的一點，n是Σ平面上rd處的單位矢量。

2 仿真方法

考慮實際情況中生物組織結(jié)構(gòu)和電導(dǎo)率相對復(fù)雜，本研究從不規(guī)則的幾何模型出發(fā)，分別建立了圓形單線圈激勵下的偏心球、正方體、橢球及電導(dǎo)率沿徑向呈余弦連續(xù)分布的球模型。在直角坐標系 x，y，z中，線圈平面平行 xoy，位于 z=0.1 m處。線圈半徑R=0.1 m，線圈截面半徑r=0.002 5 m，激勵磁場脈寬為2 μs。靜態(tài)磁場沿著 z方向均勻分布，仿真模型參數(shù)如表1所示。

表1 模型參數(shù)Tab.1 Model parameters

建立數(shù)學(xué)模型后，應(yīng)用有限元分析軟件Comsol 3.5a對模型求解域進行離散剖分，施加電磁場邊界條件，采用伽遼金加權(quán)余量有限元方法建立方程組，求解方程組得到模型中的感應(yīng)渦流。利用有限元分析結(jié)果，疊加空間均勻分布的靜態(tài)磁場，進一步得到洛倫茲力的散度即聲振動源，應(yīng)用Matlab軟件工具進行模擬聲場的計算和聲源重建工作。對聲場的仿真，假設(shè)組織的聲學(xué)特性均勻，聲波在傳播路徑上沒有散射和衰減，由式(12)求解得到模擬聲場。本研究在仿真過程中把空間固定位置處的聲壓隨著時間的變化，都相應(yīng)地轉(zhuǎn)化為聲壓隨傳播距離的變化進行計算，此處傳播距離實際為聲速與時間的乘積。同時假設(shè)聲探測器在檢測斷面上以原點為中心呈環(huán)形均勻離散分布，用來檢測不同位置接收的聲壓信號，最后應(yīng)用時間反演式(13)進行聲源圖像的重建仿真。由于實際進行成像時，很難檢測到所有位置處的聲壓信號，因此本研究最后還進行了有限角度掃描下的聲源重建仿真，分析了檢測范圍的有限性對重建聲源圖像的影響。

3 結(jié)果與分析

3.1 不同模型聲源分布

通過對不同模型的仿真，選取z=0平面顯示仿真結(jié)果。圖1為原始電導(dǎo)率在z=0平面分布，(a)～(d)分別代表偏心球、正方體、橢球、連續(xù)球模型，后續(xù)內(nèi)容中字母與圖的對應(yīng)關(guān)系與此處相同。圖2為不同模型聲源在z=0平面分布，對于電導(dǎo)率非連續(xù)分布的偏心球、正方體、橢球、連續(xù)球模型，在電導(dǎo)率躍變的邊界，感應(yīng)渦流較大，受到的洛倫茲力較大，相應(yīng)的聲振動源也大，聲源的平面分布能夠明顯地區(qū)分模型的邊界。從圖3的曲線中也可以看出，與電導(dǎo)率相同的內(nèi)部聲源相比，邊界聲源要大很多，因此，在這種強對比情況下，聲源幾乎僅能反映模型的基本幾何輪廓。對于電導(dǎo)率呈余弦分布的連續(xù)同心球模型，由于電導(dǎo)率呈緩慢地連續(xù)變化，從z=0平面聲源沿x軸的分布曲線可以看出，聲源沒有出現(xiàn)類似于前面3種模型邊界有較大跳變的情況，變化相對平緩，仿體組織邊界區(qū)分不明顯。生物體不同組織的電導(dǎo)率實際上也存在一定過渡帶，并不是嚴格躍變的，分析連續(xù)電導(dǎo)率模型對后期建立更加真實的組織模型有一定指導(dǎo)意義。

3.2 檢測聲壓隨距離的變化

圖1 原始電導(dǎo)率分布。(a)偏心球;(b)正方體;(c)橢球;(d)連續(xù)球Fig.1 Target conductivity distribution.(a)eccentric sphere;(b)cube;(c)ellipsoid;(d)continuous sphere

圖2 聲源在z=0平面分布。(a)偏心球;(b)正方體;(c)橢球;(d)連續(xù)球Fig.2 Acoustic source distribution in z=0 plane.(a)eccentric sphere;(b)cube;(c)ellipsoid;(d)continuous sphere

將位于檢測面上x軸固定點處的聲探頭接收的聲信號隨著時間的變化曲線轉(zhuǎn)化為聲壓與距離的關(guān)系，如圖4所示，橫坐標的距離代表著聲速與對應(yīng)時間的乘積。圖5為探測器接收的聲信號是以探測器位置P為原點，以R為半徑關(guān)于聲源的弧線l積分，對于電導(dǎo)率非連續(xù)分布的偏心球、正方體和橢球模型，從聲壓信號分布曲線可以看出，聲壓信號幅值發(fā)生較大變化的位置對應(yīng)著模型電導(dǎo)率變化的邊界。以偏心球為例，聲壓曲線出現(xiàn)了4個尖峰，中間的兩個尖峰之間的距離對應(yīng)著內(nèi)層球的直徑，外邊兩個尖峰的間距對應(yīng)著外層球的直徑，第一個和第二個尖峰的間距對應(yīng)著內(nèi)外層球邊界的距離。對于電導(dǎo)率呈余弦連續(xù)分布的球模型，聲源變化緩慢，因此，探測器接收的聲壓信號沒有出現(xiàn)較大的突變，變化也較為平緩。

3.3 聲源重建

圖3 聲源在z=0平面沿x軸分布。(a)偏心球;(b)正方體;(c)橢球;(d)連續(xù)球Fig.3 Acoustic source distribution along x axis in z=0 plane.(a)eccentric sphere;(b)cube;(c)ellipsoid;(d)continuous sphere

圖4 x軸正方向一個聲探頭檢測的聲信號。(a)偏心球;(b)正方體;(c)橢球;(d)連續(xù)球Fig.4 Pressure distribution with one detector located in x forward direction.(a)eccentric sphere;(b)cube;(c)ellipsoid;(d)continuous sphere

圖5 投影示意圖Fig.5 The projection diagram

由聲場中的時間反演方法即式(13)進行聲源重建得到的圖像如圖6所示。對電導(dǎo)率非連續(xù)分布模型，由于聲壓信號幅值的較大變化對應(yīng)著電導(dǎo)率變化的邊界，根據(jù)重建理論公式，反投影后邊界聲源相對內(nèi)部聲源大，使得內(nèi)外兩層球重建聲源區(qū)分不夠明顯，重構(gòu)的聲源圖像基本顯示了物體的邊界輪廓。對于電導(dǎo)率呈連續(xù)分布的球模型，聲壓信號變化平緩，反投影后的重建聲源平面圖不僅再現(xiàn)了模型的幾何形狀，同時還體現(xiàn)了聲源緩慢變化的特點。對于不連續(xù)模型，可以預(yù)先濾除幅值較大的邊界信號，通過內(nèi)外層圖像的反差來顯示邊界輪廓;而對電導(dǎo)率連續(xù)分布的模型，不存在躍變現(xiàn)象，可以通過增加采樣點獲得較好的重構(gòu)效果。由理論公式(13)可知時間反演法的實質(zhì)是投影數(shù)據(jù)沿著弧線的直接反投影重建算法，聲源圖像的像素值是所有反投影數(shù)據(jù)在該位置處的疊加，如圖5所示，聲源M和N到探測器P的距離相等，當(dāng)把聲信號反投影到非邊界聲源M時，由于聲壓信號疊加了具有較大幅值的邊界聲源N，使得反投影后M和N點得到了一個相同的投影數(shù)據(jù)疊加項，從而導(dǎo)致了重建圖像聲源的模糊效果，后期可以通過先修正后反投影的濾波反投影重建方法改善圖像質(zhì)量。

圖6 重建聲源在z=0平面分布。(a)偏心球;(b)正方體;(c)橢球;(d)連續(xù)球Fig.6 Reconstruction acoustic source distribution in z=0 plane.(a)eccentric sphere;(b)cube;(c)ellipsoid;(d)continuous sphere

3.4 探測器檢測范圍對重建圖像的影響

聲振動源較好反映了模型的幾何結(jié)構(gòu)，以重建模型基本輪廓為目的，還分析了探測器檢測范圍的有限性對重建圖像的影響。圖7是探測器以坐標原點為中心，均勻分布在圓周檢測斷面(0，π/2)范圍內(nèi)得到的重建圖像，可以看出，由于掃描角度有限，重建模型邊界有些也完全失真。分析聲源重建式(13)，得知聲壓二階導(dǎo)數(shù)前面相乘了一個比例因子，當(dāng)聲源位于探測器的法線方向時，比重最大，偏離法線方向，比重較小。因此，對于有限角度掃描的情況，能有效重建的位置取決于比例因子和探測器的掃描范圍。在環(huán)形掃描模式下，兩個測量位置的角度差為π時，投影值呈反相關(guān)系，對重建聲源是等效的，因此，在圖6中真實再現(xiàn)邊界輪廓的范圍還包括角度差為π的邊界點，該處法線方向也覆蓋了(0，π/2)范圍內(nèi)的聲探測器。

圖7 不完全數(shù)據(jù)重建聲源在z=0平面分布。(a)偏心球;(b)正方體;(c)橢球;(d)連續(xù)球Fig.7 Reconstruction acoustic source distribution in z=0 plane with incomplete data.(a)eccentric sphere;(b)cube;(c)ellipsoid;(d)continuous sphere

4 討論和結(jié)論

本研究從不規(guī)則的簡單幾何模型出發(fā)，建立了偏心球、正方體、橢球和電導(dǎo)率連續(xù)分布的球模型，應(yīng)用有限元分析方法對模型中感應(yīng)渦流進行仿真，采用Matlab軟件工具求解聲壓來模擬逆問題中的聲壓數(shù)據(jù)，通過時間反演法重建了不同模型的聲源圖像，同時還重建了有限角度掃描下的聲源圖像。從仿真結(jié)果可知，不同參數(shù)的幾何模型，探測器接收的聲信號差異較大，為獲得高質(zhì)量的重構(gòu)圖像，對于不連續(xù)模型，可以預(yù)先濾除邊界躍變來提高內(nèi)外兩層球的聲源對比度;而對電導(dǎo)率連續(xù)分布的模型，不存在躍變現(xiàn)象，可以通過增加采樣點獲得較好的重構(gòu)效果。應(yīng)用有限元和時間反演法能夠重構(gòu)聲源圖像，但本研究重建的圖像存在明顯的弧線偽跡，這與采用的直接重建算法有關(guān)，為了減少或消除偽跡對圖像質(zhì)量的影響，需要進一步尋求更高效的算法。實際情況中很難得到完備的投影數(shù)據(jù)，有限角度掃描仿真結(jié)果表明，直接采用本研究的重建算法，檢測范圍的有限性會使重建圖像失真，因此也有必要探索新的算法來完善更為接近實際成像條件下的圖像質(zhì)量。生物組織的電導(dǎo)率分布并非嚴格躍變，研究電導(dǎo)率連續(xù)分布的球模型對磁聲成像的實際應(yīng)用提供了較好的理論基礎(chǔ)，但是需要進一步建立更為真實的生物組織電導(dǎo)率連續(xù)分布模型。由于目前只是定性分析了聲源圖像的重建工作，后期對于感應(yīng)式磁聲電導(dǎo)率成像還需要深入地進行定量地探索，從而重建出能夠反映生物組織電特性的圖像。

［1］安源，任超世.電阻抗斷層功能成像技術(shù)的發(fā)展［J］.第四軍醫(yī)大學(xué)學(xué)報，2001，22(1):90－92.

［2］董秀珍.生物電阻抗成像研究的現(xiàn)狀與挑戰(zhàn)［J］.中國生物醫(yī)學(xué)工程學(xué)報，2008，27(5):641－643，649.

［3］Bayford RH.Bioimpedance tomography［J］.Annu Rev Biomed Eng，2006，8:63－91.

［4］Xu Yuan，He Bin.Magnetoacoustic tomography with magnetic induction［J］.Phys Med Biol，2005，50:5175 －5187.

［5］Li Xu，Xu Yuan，He Bin.Magnetoacoustic tomography with magnetic induction for imaging electrical impedance of biological tissue［J］.Appl Phys，2006，99:066112.

［6］Xia Rongmin，Li Xu，He Bin.Magnetoacoustic tomographic imaging of electrical impedance with magnetic induction ［J］.Applied Physics Letters，2007，91:083903.

［7］Li Xu，Xu Yuan，He Bin.Imaging Electrical impedance from acoustic measurements by means of magnetoacoustic tomography with magnetic induction(MAT-MI)［J］.IEEE Trans Biomed Eng，2007，54(2):323－330.

［8］Ma Qingyu，He Bin.Investigation on magnetoacoustic signal generation with magnetic induction and its application to electrical conductivity reconstruction ［J］. Phys Med Biol，2007，52:5085－5099.

［9］張媛.磁感應(yīng)磁聲成像的聲源產(chǎn)生機制［J］.湖州師范學(xué)院學(xué)報，2009，31(1):70－75.

［10］李珣，Li X，朱善安，等.基于時間反演方法的三維磁感應(yīng)磁聲成像電導(dǎo)率重建［J］.中國生物醫(yī)學(xué)工程學(xué)報，2009，28(1):48－52.

［11］賀文靜，劉國強，張洋，等.感應(yīng)式磁聲成像聲場正問題研究(一)—基于聲壓速度耦合方程的聲場模擬方法［J］.現(xiàn)代科學(xué)儀器，2010，(1):9－13.

［12］劉國強，賀文靜，夏慧，等.感應(yīng)式磁聲成像聲場正問題研究(二)—基于位移方程的聲場模擬方法［J］.現(xiàn)代科學(xué)儀器，2010，(1):14－16，23.

［13］劉國強，賀文靜，夏慧，等.感應(yīng)式磁聲成像聲場正問題研究(三)—基于弱形式處理洛倫茲力散度聲源的聲場模擬方法［J］.現(xiàn)代科學(xué)儀器，2010，(1):17－20.

［14］許國輝，劉志朋，李經(jīng)宇，等.感應(yīng)式磁聲成像正問題數(shù)學(xué)模型及仿真［J］.中國生物醫(yī)學(xué)工程學(xué)報，2009，28(4):481－484，489.

［15］李經(jīng)宇，殷濤，劉志朋，等.基于磁聲耦合效應(yīng)的電導(dǎo)率圖像重建研究［J］.生物醫(yī)學(xué)工程學(xué)雜志，2010，27(2):416－420.