馬虹霞 劉志朋 張順起 殷 濤

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

引言

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

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

1 理論方程

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

1.1 正問(wèn)題

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

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

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

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

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

亦即

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

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

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

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

1.2 逆問(wèn)題

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

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

2 仿真方法

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

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

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

3 結(jié)果與分析

3.1 不同模型聲源分布

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

3.2 檢測(cè)聲壓隨距離的變化

圖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

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

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

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

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

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

［2］董秀珍.生物電阻抗成像研究的現(xiàn)狀與挑戰(zhàn)［J］.中國(guó)生物醫(yī)學(xué)工程學(xué)報(bào)，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ī)制［J］.湖州師范學(xué)院學(xué)報(bào)，2009，31(1):70－75.

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

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

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

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

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

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