郭各樸 宿慧丹 丁鶴平 馬青玉

(南京師范大學(xué)物理科學(xué)與技術(shù)學(xué)院,南京 210023)

基于電阻抗層析成像的高強(qiáng)度聚焦超聲溫度監(jiān)測(cè)技術(shù)?

郭各樸 宿慧丹 丁鶴平 馬青玉?

(南京師范大學(xué)物理科學(xué)與技術(shù)學(xué)院,南京 210023)

(2017年4月8日收到;2017年6月7日收到修改稿)

作為一種對(duì)正常組織無(wú)損傷且不易引起癌細(xì)胞轉(zhuǎn)移的非入侵腫瘤治療手段,高強(qiáng)度聚焦超聲(HIFU)治療過(guò)程中焦域的溫度監(jiān)測(cè)是實(shí)現(xiàn)劑量精準(zhǔn)控制的關(guān)鍵.本文基于生物組織的溫度-電阻抗的關(guān)系,將電阻抗層析成像(EIT)和HIFU治療相結(jié)合,提出了一種利用組織焦平面的表面電壓實(shí)現(xiàn)電阻抗重構(gòu)的檢測(cè)技術(shù).建立了HIFU治療和EIT綜合系統(tǒng)模型,在考慮組織的聲吸收條件下,對(duì)三維Helmholtz方程在柱坐標(biāo)下的聲場(chǎng)計(jì)算進(jìn)行了二維簡(jiǎn)化,并引入Pennes生物熱傳導(dǎo)方程來(lái)計(jì)算HIFU焦域的聲壓和溫升分布特性;引入生物組織的溫度-電阻抗關(guān)系,基于麥克斯韋電磁場(chǎng)理論,建立了具有溫度分布HIFU焦域的電流和電壓計(jì)算模型,利用恒流注入的邊界條件實(shí)現(xiàn)電場(chǎng)計(jì)算,獲得焦平面的表面電壓分布.在數(shù)值計(jì)算中,利用實(shí)驗(yàn)聚焦換能器參數(shù),模擬了在固定聲功率下組織焦域的聲場(chǎng)和溫度場(chǎng)分布,以及中心和偏心聚焦條件下不同治療時(shí)刻的電導(dǎo)率分布;然后通過(guò)對(duì)稱(chēng)電極的循環(huán)電流注入,計(jì)算了組織模型焦平面內(nèi)的電流密度和電勢(shì)分布,獲得了焦平面圓周分布的表面電極電壓;進(jìn)一步采用修正的牛頓-拉夫遜算法,利用32×32的表面電極電壓實(shí)現(xiàn)了焦平面內(nèi)電導(dǎo)率分布的重建.結(jié)果表明,基于溫度-電阻抗關(guān)系的EIT電導(dǎo)率重建技術(shù)不但能準(zhǔn)確定位HIFU焦域中心,還能恢復(fù)HIFU治療中焦域的溫度分布,證明了EIT用于HIFU治療中溫度監(jiān)測(cè)的可行性,為其療效評(píng)估和劑量控制提供了一種無(wú)創(chuàng)電阻抗測(cè)量和成像新方法.

高強(qiáng)度聚焦超聲,電阻抗層析成像,溫度監(jiān)測(cè),表面電極電壓

1 引 言

高強(qiáng)度聚焦超聲[1?4](high intensity focused ultrasound,HIFU)是一種具有廣闊應(yīng)用前景的無(wú)創(chuàng)腫瘤治療技術(shù),它利用超聲波在組織中的穿透性和易聚焦性,將換能器發(fā)射的超聲波匯聚到腫瘤靶區(qū),利用組織的聲熱效應(yīng)產(chǎn)生65?C以上的高溫,實(shí)現(xiàn)腫瘤組織在短時(shí)間內(nèi)凝固性壞死從而達(dá)到治療腫瘤的目的.在HIFU治療過(guò)程中,既要?dú)缒[瘤細(xì)胞,又不損傷周?chē)恼＝M織,準(zhǔn)確的溫度控制是關(guān)鍵,其中實(shí)時(shí)溫度和療效監(jiān)測(cè)對(duì)HIFU的臨床應(yīng)用具有重要的意義.

為了避免在HIFU治療中插入測(cè)溫探針,減少探針對(duì)HIFU聲場(chǎng)的影響并降低癌細(xì)胞轉(zhuǎn)移概率[5],國(guó)內(nèi)外學(xué)者提出了多種無(wú)創(chuàng)測(cè)溫技術(shù).微波測(cè)溫技術(shù)利用組織溫度和熱輻射的關(guān)系,通過(guò)體外測(cè)量體內(nèi)的熱輻射來(lái)推測(cè)體內(nèi)溫度,但滲透深度有限,測(cè)量精度較差.磁共振成像測(cè)溫技術(shù)[6]通過(guò)溫度相關(guān)的擴(kuò)散系數(shù)、質(zhì)子共振頻率或弛豫時(shí)間的測(cè)量實(shí)現(xiàn)組織溫度圖像的重建,具有無(wú)創(chuàng)傷、無(wú)電離輻射、高溫度分辨率的優(yōu)點(diǎn),但其時(shí)間分辨率不高.通過(guò)將超聲測(cè)溫技術(shù)[7?11]和HIFU系統(tǒng)進(jìn)行融合,利用聲速和回波時(shí)移以及非線性等參數(shù)實(shí)現(xiàn)溫度監(jiān)測(cè),但它們的溫度變化較小,測(cè)量精度較低.近年來(lái)B超[12]被用來(lái)進(jìn)行HIFU定位引導(dǎo),監(jiān)測(cè)治療前后組織的供血變化,但是還不能實(shí)現(xiàn)高精度的溫度監(jiān)控和實(shí)時(shí)療效評(píng)價(jià).

研究表明,生物組織也是一種導(dǎo)電體,在低頻信號(hào)(＜1 MHz)激勵(lì)下正常組織的電導(dǎo)率為0—0.5 S/m[13],其導(dǎo)電能力和溫度有著明顯的對(duì)應(yīng)關(guān)系.國(guó)內(nèi)外研究發(fā)現(xiàn)生物組織的電阻抗與溫度變化存在近似線性關(guān)系[14?16],其電導(dǎo)率隨著溫度的升高而增大,其溫度-電阻抗系數(shù)(teMperature iMpedance variation factor,TIVF)約為?2%/?C,在70?C熱凝固變性時(shí)突變到?14.7%/?C.37?C和70?C時(shí)生物組織的電導(dǎo)率分別為0.41 S/m和0.79 S/m,其變化幾乎達(dá)到100%.和聲阻抗相比,生物組織的電阻抗范圍和變化率更高,這為無(wú)創(chuàng)測(cè)溫提供了物理基礎(chǔ).

在本實(shí)驗(yàn)室近期的研究中,利用組織模型的電阻抗相對(duì)變化[17](relative iMpedance variation,RIV)進(jìn)行了HIFU焦域的電阻抗測(cè)量,結(jié)果表明當(dāng)HIFU聲功率一定時(shí),組織模型的RIV和治療時(shí)間呈線性關(guān)系;在達(dá)到HIFU治療療效時(shí)(焦域徑向±0.4 mm區(qū)域達(dá)到70?C),組織模型的RIV和HIFU聲功率呈現(xiàn)反比關(guān)系,而所需治療時(shí)間和聲功率的平方成反比,證明組織模型的RIV可以用來(lái)實(shí)現(xiàn)HIFU治療過(guò)程中組織焦域的溫度監(jiān)測(cè)和療效評(píng)估.但是由于HIFU焦域的尺寸很小,其電阻抗變化對(duì)模型RIV的影響較小,因此對(duì)電阻抗測(cè)量系統(tǒng)的靈敏度提出了更高的要求.同時(shí),由于RIV測(cè)量的是HIFU治療過(guò)程中組織模型的整體電阻抗變化,很難實(shí)現(xiàn)高精度的焦點(diǎn)定位和焦域溫度反演,進(jìn)一步將HIFU治療過(guò)程中的溫度-電阻抗關(guān)系和焦域的準(zhǔn)確定位相結(jié)合,進(jìn)行電阻抗/溫度圖像的重建在HIFU治療的溫度監(jiān)測(cè)中具有重大的研究?jī)r(jià)值和應(yīng)用前景.

近年來(lái),電阻抗成像技術(shù)(electrical iMpedance tomography,EIT)[18?20]得到了迅速的發(fā)展,并在生物醫(yī)學(xué)領(lǐng)域得到了廣泛的應(yīng)用.EIT根據(jù)生物組織內(nèi)部電阻抗分布的不同,通過(guò)對(duì)物體表面電流和電壓的測(cè)量來(lái)重建物體內(nèi)部電導(dǎo)率的分布.由于EIT技術(shù)不使用射線,激勵(lì)電流在安全范圍以內(nèi),對(duì)人體無(wú)害,成本低廉,且可以實(shí)現(xiàn)HIFU治療聲學(xué)系統(tǒng)和電阻抗測(cè)量的電學(xué)系統(tǒng)的完全分離,減少相互干擾,提高測(cè)量的可靠性,因此用EIT來(lái)測(cè)量HIFU焦域的電阻抗分布是一種新方法.

本文基于組織的溫度-電阻抗關(guān)系,建立了HIFU治療和EIT測(cè)量系統(tǒng)模型,研究了HIFU的聲場(chǎng)和溫度場(chǎng)分布特性,分析了HIFU焦域的電阻抗分布,提出了一種基于EIT的HIFU焦域溫度監(jiān)測(cè)技術(shù).首先,利用有限元建立了系統(tǒng)的仿真模型,計(jì)算了HIFU過(guò)程中組織模型的聲場(chǎng)和溫度場(chǎng)分布;然后基于組織的溫度-電阻抗關(guān)系,將聲學(xué)系統(tǒng)擴(kuò)展到電學(xué)系統(tǒng),利用焦平面表面對(duì)稱(chēng)電極的循環(huán)電流注入,計(jì)算HIFU焦域的電流密度和電勢(shì)分布,得到焦平面上的表面電極電壓;最后采用修正的牛頓-拉夫遜算法,利用模擬的32×32表面電極電壓進(jìn)行焦平面的電導(dǎo)率分布重建.中心和偏心聚焦條件下不同治療時(shí)間的仿真結(jié)果表明,重建的電導(dǎo)率分布能較好地反映出HIFU的焦點(diǎn)位置和焦域的電導(dǎo)率變化,間接反映了HIFU治療過(guò)程中焦域的溫升情況.研究結(jié)果證明了EIT用于HIFU焦域電阻抗分布重建的可行性,為HIFU治療中焦域的精確定位和溫度監(jiān)測(cè)以及療效評(píng)估提供了一種電阻抗測(cè)量和成像的新方法.

2 原理與方法

Crum等[21]對(duì)超聲非線性對(duì)HIFU焦域溫升的研究表明,由于熱傳導(dǎo)的影響,超聲非線性對(duì)HIFU焦域的組織溫升影響較小;另外Myers和Soneson[22?24]證明n階諧波所產(chǎn)生的溫升和其諧波熱源的比值為log(n)/n,在HIFU治療中非線性效應(yīng)產(chǎn)生明顯焦點(diǎn)溫升誤差(達(dá)到10%)的聲功率臨界值為115W.本研究所用的聲功率為15.68W時(shí),非線性溫升誤差遠(yuǎn)小于2%,可以采用線性模型進(jìn)行聲場(chǎng)計(jì)算和溫升估計(jì),為進(jìn)一步HIFU焦域的電壓測(cè)量和EIT重建提供基礎(chǔ).

考慮到媒質(zhì)的黏滯性對(duì)聲波的吸收,聚焦超聲換能器所激發(fā)的聲波在生物組織內(nèi)部傳播時(shí)會(huì)產(chǎn)生能量衰減和溫度升高,其聲壓p滿足一維波動(dòng)方程[25?27]:

圖1 HIFU治療和EIT測(cè)量系統(tǒng)原理圖Fig.1.Sketch Map of the coMp rehensive systeMof H IFU therapy and EIT Measu reMent.

HIFU治療系統(tǒng)的原理如圖1(a)所示,聚焦換能器表面任意超聲陣元的波動(dòng)方程可以簡(jiǎn)化為Helmholtz[26]方程:

由于換能器聲場(chǎng)的軸對(duì)稱(chēng)性,(3)式可以在二維軸對(duì)稱(chēng)圓柱坐標(biāo)下簡(jiǎn)化為

其中,er和ez分別是沿r和z方向的單位矢量.將cc= ω/k′代入(4)式,將公式左右兩邊同乘(?r),得到黏滯介質(zhì)中的簡(jiǎn)化方程:

眾所周知,超聲的聲熱效應(yīng)[26?28]是超聲熱療的基本原理.一定強(qiáng)度的超聲在組織媒介中傳播時(shí),部分聲能被組織吸收轉(zhuǎn)化為熱能,這種被吸收的聲能就是HIFU治療的熱源.假設(shè)聲波在沿z軸傳播過(guò)程中通過(guò)一個(gè)面積為d S、長(zhǎng)度為d z的體積元媒質(zhì),如z處的聲強(qiáng)為I1,z+d z處的聲強(qiáng)為I2,考慮到媒質(zhì)對(duì)聲波能量的吸收,可知I1＞I2.在一維聲傳播條件下,聲場(chǎng)中單位時(shí)間內(nèi)單位體積的媒質(zhì)所吸收的聲能量可以由聲強(qiáng)表示:

在三維空間的聲傳播中,熱源則可以由聲強(qiáng)的空間梯度[28]計(jì)算,

在平面波近似下,聲波沿z方向傳播,若入射超聲強(qiáng)度I0,傳播距離z處的聲強(qiáng)可以表示為I=I0exp(?2αaz),其中αa為吸收系數(shù),代入(7)式得到[28]:

在生物組織中,入射聲能量的損失一般由聲衰減系數(shù)α表征,其為吸收系數(shù)αa和散射系數(shù)αs的總和.實(shí)際計(jì)算中很難區(qū)分組織對(duì)聲能量的吸收和散射分量,一般直接將衰減系數(shù)等效為吸收系數(shù)[28],即α≈αa.在忽略介質(zhì)中的熱耗散和熱傳導(dǎo),傳播距離z處的超聲熱源表示為[28]

在聲場(chǎng)計(jì)算中,HIFU焦點(diǎn)處的聲強(qiáng)可以由聲壓各次諧波幅值表示[28]:

其中,ρt和ct分別為組織密度和聲速,Cn是n次諧波的聲壓幅值.在本研究中所用超聲功率不大,忽略非線性影響,(10)式可以簡(jiǎn)化為用基波聲壓計(jì)算.

在不考慮血液流動(dòng)影響的前提下,將組織所吸收的熱量應(yīng)用到Pennes生物熱傳導(dǎo)方程[28,29]中,

其中Kt是組織熱導(dǎo)率,T為組織溫度,T0為初始溫度.由于HIFU的熱效應(yīng),在焦點(diǎn)處形成中心溫度最高、周?chē)鷾囟容^低、具有明顯的溫度梯度分布的橢球狀(mm)的溫升焦域.為了定量分析組織電阻抗隨溫度的變化,將TIVF[14?16]應(yīng)用到電導(dǎo)率中,得到組織的溫度-電導(dǎo)率分段函數(shù)為

可見(jiàn),隨著HIFU治療時(shí)間的增長(zhǎng),生物組織焦域的溫度逐漸升高,電導(dǎo)率隨溫度的升高而增大,在70?C組織凝固時(shí)產(chǎn)生快速的電導(dǎo)率提升,在焦域內(nèi)形成電導(dǎo)率的梯度分布,這為EIT在HIFU中的溫度監(jiān)測(cè)提供了物理基礎(chǔ).

基于EIT的HIFU焦域溫度監(jiān)測(cè)技術(shù)的正問(wèn)題研究包括聲場(chǎng)、溫度和電導(dǎo)率分布,焦平面內(nèi)的電勢(shì)分布以及邊界電極的電壓分布.圖1(b)顯示了HIFU治療中組織模型(半徑為R)焦平面的電導(dǎo)率分布,在邊界圓周上均勻設(shè)置N個(gè)電極,用對(duì)稱(chēng)電極1/17電流注入,計(jì)算焦平面內(nèi)的電流密度和電勢(shì)分布,得到32個(gè)邊界電極的電壓值;然后利用對(duì)稱(chēng)電極的循環(huán)電流注入,重復(fù)以上過(guò)程,得到邊界電極電壓Vij(i,j=1,2,...,N),其中i和j表示激勵(lì)電極和測(cè)量電極的序號(hào).

EIT正問(wèn)題求解是在已知HIFU焦域的內(nèi)部電阻抗分布的前提下,且具有特殊邊界條件的電場(chǎng)計(jì)算[30,31],組織模型可以等效為導(dǎo)體,其電場(chǎng)分布滿足麥克斯韋方程組[32],在10—100 kHz的低頻電流激勵(lì)下,忽略介電常數(shù)的影響,得到場(chǎng)域的表達(dá)式[33]:

其中,σ(T)為場(chǎng)域內(nèi)電導(dǎo)率分布,T為節(jié)點(diǎn)的溫度;Γ1為第一類(lèi)邊值條件,表示已知位函數(shù)在場(chǎng)域邊界上各點(diǎn)的值,φ為場(chǎng)域內(nèi)電勢(shì)分布函數(shù),φ0為給定的邊界電位,初始時(shí)認(rèn)為給定點(diǎn)邊界電位為0;Γ2為第二類(lèi)邊值問(wèn)題,表示已知位函數(shù)在場(chǎng)域邊界上各點(diǎn)的法向?qū)?shù)值;Jn是給定邊界注入的電流密度.對(duì)(13)式采用變分法求解[33],建立拉普拉斯方程的泛函

根據(jù)格林定理,考慮模型的邊界條件,并使整個(gè)場(chǎng)域的總能量最小,得到

在HIFU治療中,焦域的電導(dǎo)率分布隨著溫升而改變,結(jié)合激勵(lì)電極的邊界條件,利用有限元算法可以計(jì)算出組織模型內(nèi)各剖分單元和節(jié)點(diǎn)的電壓,進(jìn)一步獲得邊界電極測(cè)量電壓Vij.

基于EIT的HIFU焦域溫度監(jiān)測(cè)的逆問(wèn)題是通過(guò)已經(jīng)獲得的模型邊界電極電壓和電流激勵(lì)模式,利用Vij重建出模型內(nèi)的電阻抗分布.修正的牛頓-拉夫遜(MNR)算法[18,34,35]是通過(guò)不斷迭代來(lái)改變電阻率分布,進(jìn)而使目標(biāo)函數(shù)(重建的邊界電極電壓和測(cè)量電壓之間誤差范數(shù)的平方)最小.如組織的電導(dǎo)率為σ(T),則相應(yīng)的電阻率分布為T(mén).在HIFU治療的正問(wèn)題中,通過(guò)對(duì)稱(chēng)電極的循環(huán)電流注入得到32×32邊界電極測(cè)量電壓Vij;在MNR電阻抗重建的逆問(wèn)題中,假設(shè)相應(yīng)對(duì)稱(chēng)電極電流注入時(shí)計(jì)算得到邊界電極電壓為Uij(ρ),則目標(biāo)函數(shù)為

其中N=32.通過(guò)不斷迭代使目標(biāo)函數(shù)最小,得到焦平面的穩(wěn)定電阻率分布.為使目標(biāo)函數(shù)最小,令f′(ρ)=U′(ρ)(U(ρ)?V)=0,其中U′(ρ)稱(chēng)為Jacobian矩陣[36].對(duì)f′(ρ)泰勒級(jí)數(shù)展開(kāi),并保留線性項(xiàng)得到f′≈ f′(ρk)+f′′(ρk)?ρk,其中?ρk=ρk+1? ρk,f′′(ρk)=[J(ρk)]TJ(ρk)為Hessian矩陣,?ρk= ?[[J(ρk)]TJ(ρk)]?1[J(ρk)]T[U(ρk)?V],得到MNR的迭代公式為

其中k是迭代次數(shù). 為了避免對(duì)病態(tài)[J(ρk)]TJ(ρk)求逆,在重建中引入Tikhonov正則化[33]修正,并將其補(bǔ)償項(xiàng)表示為ρ的平方函數(shù)形式,得到重構(gòu)電阻率分布的迭代公式為

將MNR重建方法應(yīng)用到HIFU焦域的電阻率重建中,先基于EIT正問(wèn)題得到的邊界電極測(cè)量電壓Vij,假設(shè)模型焦平面內(nèi)電阻率分布均勻,利用有限元模型計(jì)算得到邊界電極電壓Uij(ρ),建立目標(biāo)函數(shù)f(ρ),通過(guò)求解不同電阻率分布下的雅克比矩陣和海森矩陣,獲得新的迭代電阻率分布,直到目標(biāo)函數(shù)小于預(yù)設(shè)值,此時(shí)的電阻率分布ρk即為重建出的組織模型焦平面的電阻率分布,進(jìn)一步利用σ(T)=1/ρk可以計(jì)算出電導(dǎo)率分布.

3 數(shù)值計(jì)算

如圖1所示,由于換能器和聲傳播的對(duì)稱(chēng)性,在有限元[37]計(jì)算中,HIFU聲場(chǎng)和溫度場(chǎng)以及組織電導(dǎo)率分布采用二維軸對(duì)稱(chēng)的柱坐標(biāo)模型,其中軸向z是聲傳播方向,r是半徑方向,仿真區(qū)域?yàn)槌晸Q能器、水域環(huán)境以及組織模型.通過(guò)調(diào)整換能器表面振速控制輸出聲功率,計(jì)算組織焦域的聲場(chǎng)、溫度場(chǎng)以及電導(dǎo)率分布.為了在保證計(jì)算精度的前提下提高計(jì)算速度,水域環(huán)境區(qū)域剖分網(wǎng)格尺寸為λWater/4,組織區(qū)域和焦域的剖分網(wǎng)格尺寸分別λTissue/4和λTissue/8.設(shè)置聚焦超聲換能器的直徑和焦距均為10 cm,中心頻率1.13 MHz.直徑和高度分別為32mm和35 mm的圓柱形組織模型中心放在超聲換能器的焦域處.在模型焦平面的表面均勻設(shè)置32個(gè)電極測(cè)量表面電極的電壓.在HIFU治療中,隨著超聲的作用,模型內(nèi)焦域的溫度升高,其電導(dǎo)率隨之提高,形成電導(dǎo)率的梯度分布.模型采用仿組織透明凝膠[38],其物理參數(shù)和人體組織較為接近.計(jì)算中水和凝膠組織以及人體組織的相關(guān)參數(shù)列于表1.結(jié)合HIFU系統(tǒng)的實(shí)驗(yàn)參數(shù),將換能器表面振幅設(shè)定為10.2 nm,通過(guò)有限元計(jì)算得到焦域處的聲壓以及聲強(qiáng)分布,定義焦平面內(nèi)聲壓衰減6 dB的面積為有效截面面積,通過(guò)I(x,y)d A計(jì)算得到聲功率為15.68 W[17],并以此作為聲源參數(shù)進(jìn)行相關(guān)計(jì)算.

表1 溫度為293 K(20?C)時(shí)的仿真參數(shù)Tab le 1.ParaMeters used in siMu lation at the temperature of 293 K(20?C).

通過(guò)有限元仿真得到固定聲功率15.68 W時(shí)HIFU焦域的軸向剖面和徑向焦平面的聲壓分布,結(jié)果如圖2(a)和圖2(b)所示,可見(jiàn)HIFU焦域呈現(xiàn)橢球狀,焦點(diǎn)處的聲壓最大,隨著軸向和徑向范圍的擴(kuò)大,聲壓逐漸降低.結(jié)合(11)式,計(jì)算不同治療時(shí)間(?t)組織焦域的溫升,得到如圖3所示的焦平面溫度分布.在?t=1 s時(shí),焦域半徑小于0.5 mm,焦域中心的溫度較低,未達(dá)到70?C.隨著治療時(shí)間的延長(zhǎng),焦域的能量逐漸積累,焦點(diǎn)及周?chē)M織的溫度不斷升高,同時(shí)由于組織的熱擴(kuò)散,焦域面積不斷擴(kuò)大,在焦平面上形成圓形的焦斑.在?t=2 s時(shí),焦域中心0.2 mm半徑范圍內(nèi)的溫度超過(guò)70?C,達(dá)到了治療效果.在?t=3 s時(shí),焦域半徑超過(guò)1.3 mm,其中超過(guò)在半徑0.5 mm的范圍內(nèi)溫度超過(guò)70?C.進(jìn)一步將組織的溫度-電阻抗關(guān)系應(yīng)用到圖3(a)中,得到如圖3(b)所示的不同治療時(shí)間的焦平面電導(dǎo)率分布.可見(jiàn)隨著的延長(zhǎng),HIFU焦域的溫度逐漸升高,電導(dǎo)率相應(yīng)提高,形成焦域中心導(dǎo)電能力強(qiáng),周?chē)鷮?dǎo)電性能弱的梯度分布.

圖2 (網(wǎng)刊彩色)固定聲功率為15.68W時(shí),HIFU的(a)軸向剖面和(b)焦平面的聲壓分布Fig.2.(color on line)Pressu re distributions of(a)axial p rofi le and(b)focal p lane for HIFU at a fi xed acoustic power of 15.68 W.

圖3 (網(wǎng)刊彩色)不同治療時(shí)刻HIFU焦平面的(a)溫度和(b)電導(dǎo)率分布及(c)偏心聚焦的電導(dǎo)率分布Fig.3.(color on line)Focal d istribu tions of(a)teMperatu re and(b)conductivity for centric HIFU therapy,as well as(c)the conductivity distributions for eccentric HIFU therapy at diff erent treatMent tiMes.

為了證明研究的普適性,將焦點(diǎn)沿著徑向移動(dòng)8 mm形成HIFU偏心聚焦,得到了如圖3(c)所示的不同治療時(shí)間HIFU焦平面的電導(dǎo)率分布.由于組織媒質(zhì)相同,模型內(nèi)不同位置的聲學(xué)特性和電阻抗特性相同,因此除了焦域位置的移動(dòng),焦平面上焦點(diǎn)處的電導(dǎo)率分布和圖3(b)基本一致.

在HIFU治療的同時(shí),通過(guò)對(duì)稱(chēng)電極的電流注入,在模型內(nèi)部形成相應(yīng)的電場(chǎng)分布,其變化能反映焦域的溫度變化.圖4(a)和圖4(b)分別顯示了HIFU中心和偏心聚焦時(shí),在電極1/17電流注入,?t=1,2,3 s時(shí)焦平面內(nèi)的電流密度和電勢(shì)分布.如圖4(a)的虛線環(huán)所示,隨著治療時(shí)間的延長(zhǎng),焦域的溫度和電導(dǎo)率升高,原來(lái)較為均勻分布的電流向焦域的中心匯聚,形成向上彎曲的電勢(shì)等位線.對(duì)于如圖4(b)所示的偏心聚焦,在焦域內(nèi)也形成匯聚的電流和彎曲的電勢(shì)等位線,其彎曲程度由焦域的電導(dǎo)率分布決定.

圖4 (網(wǎng)刊彩色)對(duì)稱(chēng)電極1/17電流注入,HIFU(a)中心和(b)偏心聚焦時(shí)焦平面的電流密度和電勢(shì)分布Fig.4.(color on line)D istributions of cu rrent density and electrical potential of the focal p lane With(a)centric and(b)eccentric focusing for the current in jection froMthe electrodes 1/17.

為了進(jìn)行基于EIT的電阻抗分布重建,采用如圖1所示的三角形網(wǎng)格剖分[39],其中剖分節(jié)點(diǎn)數(shù)為1225,剖分單元數(shù)為2320,并將有限元計(jì)算得到的焦平面電導(dǎo)率分布和和電勢(shì)分布導(dǎo)入到所建立的剖分網(wǎng)格中,將組織模型的焦平面圓周32等分,獲得循環(huán)對(duì)稱(chēng)電極電流注入時(shí)的邊界電極測(cè)量電壓Vij.圖5(a)和圖5(b)分別顯示了對(duì)稱(chēng)電極1/17和9/25電流注入時(shí)32個(gè)邊界電極上的測(cè)量電壓分布.可見(jiàn)隨著治療時(shí)間的延長(zhǎng),焦域溫度升高,電阻率降低,相同電流注入時(shí)邊界電極電壓降低.在改變電流注入電極后,邊界電極的電壓分布發(fā)生相應(yīng)的變化.然而,由于相對(duì)于模型尺寸,焦域面積很小,電阻抗差異較小,因此電極電壓分布隨著治療時(shí)間和溫度以及位置的變化不明顯.圖5(a)和圖5(b)的內(nèi)插局部放大圖分別顯示了兩種情況下測(cè)量得到電極7/8和17/18的電壓分布,可以看出相應(yīng)電極的電壓產(chǎn)生了幅度較小的變化(mV量級(jí)).

圖5 (網(wǎng)刊彩色)對(duì)稱(chēng)電極(a)1/17和(b)9/25電流注入時(shí),模型焦平面表面電極的電壓分布Fig.5.(color on line)E lectrical voltage distributions of the surface electrodes in the focal p lane With the current in jections froMthe symMetric electrodes of 1/17 and 9/25.

圖6 (網(wǎng)刊彩色)H IFU中心聚焦時(shí),(a)模型焦平面的網(wǎng)格化電導(dǎo)率分布和(b)重建的電導(dǎo)率分布Fig.6.(color on line)(a)G rid d istributions and(b)reconstructed distributions of conductivity in the focal p lane for centric focusing.

在EIT的電阻抗重建中,根據(jù)HIFU治療中組織焦域的電導(dǎo)率改變,設(shè)定初始均勻電導(dǎo)率為0.73 S/m,利用Matlab編程計(jì)算焦平面上的電極電位分布Uij(ρ),建立目標(biāo)函數(shù),并通過(guò)求解每次迭代電阻率分布下的雅克比矩陣和海森矩陣獲得新的迭代電阻率分布,當(dāng)?shù)螖?shù)達(dá)到30次時(shí),獲得了較穩(wěn)定的重建結(jié)果.圖6(a)和圖6(b)分別顯示了HIFU中心聚焦時(shí)模型焦平面的網(wǎng)格化電導(dǎo)率分布和基于邊界電極電壓的重建電導(dǎo)率分布,為了進(jìn)行重建效果比較,HIFU偏心聚焦時(shí)的相應(yīng)結(jié)果在圖7(a)和圖7(b)中給出.可見(jiàn)重建結(jié)果能夠較為準(zhǔn)確地反映出治療區(qū)域的位置和形狀以及電導(dǎo)率的變化趨勢(shì),能夠?qū)崿F(xiàn)準(zhǔn)確的焦點(diǎn)定位,重建焦域的尺寸和電導(dǎo)率隨治療時(shí)間的增長(zhǎng)而變大,和模型計(jì)算結(jié)果的變化趨勢(shì)一致,能夠反映出不同治療時(shí)間HIFU焦域的電導(dǎo)率差異.但是由于中心聚焦的焦域具有完全對(duì)稱(chēng)性,其邊界電極電壓的偏差較小,因此重建焦域的精度不如偏心聚焦效果好.

由于在電導(dǎo)率的網(wǎng)格化處理中使用了三角形內(nèi)部的數(shù)據(jù)平均,所得到的最大電導(dǎo)率小于有限元的仿真結(jié)果.同時(shí)由于重建算法的病態(tài)性,所得到的電導(dǎo)率小于實(shí)際值.為了比較重建效果,將重建結(jié)果用對(duì)應(yīng)時(shí)刻模型電導(dǎo)率的最大值進(jìn)行歸一化處理,得到如圖8所示的不同治療時(shí)間模型焦平面的網(wǎng)格化電導(dǎo)率徑向分布和重建電導(dǎo)率的徑向分布.雖然重建結(jié)果比實(shí)際電導(dǎo)率范圍寬,但是仍然可以精確定位焦點(diǎn)位置,反映焦平面內(nèi)焦域溫度和電導(dǎo)率的分布;隨著治療時(shí)間的增加,焦域的電導(dǎo)率變化增大,重建的電導(dǎo)率分布和模型電導(dǎo)率分布趨勢(shì)一致,能夠較好地反映組織電導(dǎo)率隨溫度的變化特性.

圖7 (網(wǎng)刊彩色)HIFU偏心聚焦時(shí),(a)模型焦平面的網(wǎng)格化電導(dǎo)率分布和(b)重建的電導(dǎo)率分布Fig.7.(color on line)(a)G rid d istributions and(b)reconstructed distributions of conductivity in the focal p lane for eccentric focusing.

圖8 (網(wǎng)刊彩色)H IFU(a)中心和(b)偏心聚焦時(shí),不同治療時(shí)間組織模型焦平面的網(wǎng)格化電導(dǎo)率歸一化徑向分布和重建電導(dǎo)率的歸一化徑向分布Fig.8.(color on line)NorMalized rad ial grid d istributions and reconstructed d istribu tions of conductivity in the focal p lane for(a)centric and(b)eccentric focusing.

4 討 論

本研究中由于HIFU的聲功率不太高,超聲非線性所產(chǎn)生溫升對(duì)線性聲場(chǎng)所產(chǎn)生的誤差較小,為了簡(jiǎn)化溫度場(chǎng)的計(jì)算,選用黏滯媒質(zhì)中的線性方程求解,取得較好的計(jì)算效果.但實(shí)際中,HIFU聲場(chǎng)一般為非線性,要獲得更為準(zhǔn)確的聲場(chǎng)和溫度場(chǎng),甚至考慮聲空化的影響,在進(jìn)一步的研究中采用Westervelt[40],KZK[41,42]和SBE[43]方程求解,會(huì)得到更為精確HIFU焦域的溫度分布.

理論和模擬結(jié)果證明,雖然測(cè)量的邊界電極的電壓變化較小,但利用EIT可以精準(zhǔn)定位HIFU焦域的位置和焦域內(nèi)電阻抗的變化趨勢(shì),重建圖像中電導(dǎo)率的差異可以反映焦域的溫度分布特性,在HIFU的療效監(jiān)測(cè)中有著廣泛的應(yīng)用前景.另外,當(dāng)前EIT逆問(wèn)題求解算法主要包括動(dòng)態(tài)重構(gòu)算法(例如等位線反投影法)和靜態(tài)重構(gòu)算法(如牛頓-拉夫遜迭代算法)等.鑒于動(dòng)態(tài)重構(gòu)算法計(jì)算速度快和效率高的優(yōu)點(diǎn),本研究首先采用等位線投影法對(duì)32電極的測(cè)量數(shù)據(jù)進(jìn)行了仿真計(jì)算,然而由于HIFU焦域的面積小,電導(dǎo)率變化范圍大,而邊界電極的電壓變化小,因此重建算法對(duì)電阻抗的變化不敏感,重建圖像的精度差.雖然修正的牛頓-拉夫遜算法是一種靜態(tài)的重建算法,但具有良好的收斂性,重建結(jié)果較動(dòng)態(tài)算法更好,但由于該算法重建問(wèn)題的病態(tài)性會(huì)隨剖分網(wǎng)格的增加而增加,電阻抗重建的網(wǎng)格剖分不能過(guò)密,限制了重建精度的提高,因此需要進(jìn)一步研究EIT重建算法,在較少的網(wǎng)格剖分和表面電極條件下,獲得更高精度的重建圖像.

本研究使用循環(huán)對(duì)稱(chēng)電極電流注入,其電流密度和電勢(shì)分布對(duì)位置較為敏感,當(dāng)HIFU焦域在模型中心時(shí),不同位置激勵(lì)下的表面電極電壓的差異較小,因此電阻抗重建的精度和分辨率較低.為了獲得較大的表面電極電壓差異,需要改進(jìn)電極注入方式,實(shí)現(xiàn)了多種激勵(lì)模式下的電阻抗分布重建.當(dāng)前EIT理論較為成熟,算法簡(jiǎn)便,設(shè)備廉價(jià),對(duì)生物組織傷害小,重建速度較快,本研究把EIT和HIFU相結(jié)合,實(shí)現(xiàn)電學(xué)系統(tǒng)和聲學(xué)系統(tǒng)完全分離,對(duì)HIFU治療系統(tǒng)干擾小,可以同步實(shí)現(xiàn)HIFU治療和電阻抗測(cè)量,利用組織溫度-電導(dǎo)率的關(guān)系,可以實(shí)現(xiàn)HIFU治療過(guò)程中焦域溫度的反演.進(jìn)一步研究中,將模型對(duì)稱(chēng)電極的RIV測(cè)量和邊界電極電壓監(jiān)測(cè)相結(jié)合,可以有效提高焦點(diǎn)定位和焦域溫度監(jiān)測(cè)的準(zhǔn)確性,為HIFU療效評(píng)估和劑量控制提供實(shí)時(shí)電阻抗測(cè)量和成像新方法.

5 結(jié) 論

本文針對(duì)HIFU治療中的無(wú)創(chuàng)溫度監(jiān)測(cè)難題,基于組織的溫度-電阻抗關(guān)系,提出了一種基于EIT的HIFU焦域的溫度和療效監(jiān)測(cè)技術(shù).本研究將HIFU治療的聲學(xué)系統(tǒng)和電阻抗測(cè)量的電學(xué)系統(tǒng)有機(jī)結(jié)合,利用組織模型表面的電阻抗測(cè)量來(lái)實(shí)現(xiàn)HIFU焦域的電阻抗分布的重建,實(shí)現(xiàn)溫度和療效估計(jì).建立了HIFU治療和EIT測(cè)量系統(tǒng)模型,模擬了固定聲功率(15.68W)時(shí)組織焦域的聲場(chǎng)和溫度場(chǎng)分布,以及中心和偏心聚焦的電導(dǎo)率分布.通過(guò)對(duì)稱(chēng)電極的旋轉(zhuǎn)電流注入,計(jì)算了圓柱組織模型焦平面內(nèi)的電流密度和電勢(shì)分布,獲得了焦平面的表面電極電壓分布.在電阻抗重建中,采用修正的牛頓-拉夫遜算法,通過(guò)迭代計(jì)算實(shí)現(xiàn)了焦平面的電導(dǎo)率分布重建.結(jié)果表明,重建的電導(dǎo)率分布能準(zhǔn)確反映HIFU焦域的位置和范圍,體現(xiàn)隨治療時(shí)間增長(zhǎng)的電導(dǎo)率變化,證明EIT在溫度監(jiān)測(cè)中應(yīng)用的可行性.本研究為HIFU治療中焦點(diǎn)的精確定位,焦域的實(shí)時(shí)溫度監(jiān)測(cè)和療效評(píng)估以及劑量控制提供了一種無(wú)創(chuàng)電阻抗測(cè)量和成像新方法.

[1]Hutchinson L 2011 Nat.Rev.C lin.Oncol.8 385

[2]K ennedy J E 2005 Nat.Rev.Canc.5 321

[3]Q ian S Y,Wang H Z 2001 Acta Phys.Sin.50 501(in Chinese)[錢(qián)盛友,王鴻樟 2001物理學(xué)報(bào) 50 501]

[4]Gavrilov L R 2013 J.Acoust.Soc.AMer.133 4348

[5]Jiang L X,Hu B 2006 Tech.Acoust.25 43(in Chinese)[姜立新,胡兵2006聲學(xué)技術(shù) 25 43]

[6]Shen J,Shen J L,Zou J Z 2007 J.U ltrasound C lin.Med.9 486(in Chinese)[沈潔,申俊玲,鄒建中2007臨床超聲醫(yī)學(xué)雜志9 486]

[7]Ye G,SMith P P,Noble J A 2010 U ltrasound Med.Biol.36 234

[8]Daniels MJ,Varghese T,Madsen E L,Zagzebski J A 2007 Phys.Med.Biol.52 4827

[9]Fan T B,Zhang D,Zhang Z,Ma Y,Gong X F 2008 Chin.Phys.B 17 3372

[10]Anand A,KaczkoWski P J 2004 J.Acoust.Soc.AMer.115 2490

[11]Ma Y,Zhang D,Gong X F,Liu X Z,Ma Q Y,Q iu Y Y 2007 Chin.Phys.16 2745

[12]Fan L Z,Luo F,Yu D Y,Liu X T,Zhang J,X ie MX 2005 C lin.J.Med.Instru.29 115(in Chinese)[范良志,羅飛,喻道遠(yuǎn),劉夏天,張靜,謝明星 2005中國(guó)醫(yī)療器械雜志29 115]

[13]Gabriel C,PeyMan A,G rant E H 2009 Phys.Med.Biol.54 4863

[14]Zurbuchen U,HolMer C,LehMann K S,Stein T,Roggan A,Seifarth C,Buh r H J,Ritz J P 2010 In t.J.HypertherMia 26 26

[15]G riffi ths H,AhMed A 1987 C lin.Phys.Physiol.Meas.8 147

[16]Cai H,You F S,Shi X T,Fu F,Liu R G,Tang C,Dong X Z 2010 Chin.Med.Equip.J.31 8(in Chinese)[蔡華,尤富生,史學(xué)濤,付峰,劉銳崗,湯池,董秀珍2010醫(yī)療衛(wèi)生裝備31 8]

[17]Su H D,Guo G P,Ma Q Y,Tu J,Zhang D 2017 Chin.Phys.B 26 054302

[18]Li K Q 2015 M.S.Dissertation(Nan jing:Nan jing University of Posts and TelecomMunications)(in Chinese)[李凱強(qiáng)2015碩士學(xué)位論文 (南京:南京郵電大學(xué))]

[19]Xu G X 2004 Ph.D.D issertation (Chongqing:Chongqing University)(in Chinese)[徐管鑫 2004博士學(xué)位論文(重慶:重慶大學(xué))]

[20]Zhang L 2013 Ph.D.D issertation(Nan jing:Nan jing University of Science and Technology)(in Chinese)[張麗2013博士學(xué)位論文(南京:南京理工大學(xué))]

[21]Curra F P,Mourad P D,Khokh lova V A,CruML A 2000 IEEE Trans.U ltrason.Ferroelect.Freq.Con trol 47 1077

[22]Soneson J E,Myers MR 2010 IEEE Trans.U ltrason.Ferroelect.Freq.Con trol 57 2450

[23]Myers MR,Soneson J E 2009 J.Acoust.Soc.AMer.126 425

[24]Soneson J E,Myers MR 2007 J.Acoust.Soc.AMer.122 2526

[25]Du G H,Zhu Z M,Gong X F 2012 FundaMentals of Acoustics(Nan jing:Nan jing University Press)pp292–305(in Chinese)[杜功煥,朱哲民,龔秀芬 2012聲學(xué)基礎(chǔ)(南京:南京大學(xué)出版社)第292—305頁(yè)]

[26]Cheng J C 2013 Principles of Acoustics(Beijing:Science Press)pp571–576(in Chinese)[程建春 2013聲學(xué)原理 (北京:科學(xué)出版社)第571—576頁(yè)]

[27]Wan MX,Zong J Y,Wang S P 2010 BioMedical U ltrasound(Beijing:Science Press)pp649–669(in Chinese)[萬(wàn)明習(xí),宗瑜瑾,王素品2010生物醫(yī)學(xué)超聲學(xué)(北京:科學(xué)出版社)第649—669頁(yè)]

[28]Zhang D,Guo X S,Ma Q Y,Tu J 2014 FundaMen ta ls of Medical U ltrasound(Beijing:Science Press)pp415–418(in Chinese)[章東,郭霞生,馬青玉,屠娟 2014醫(yī)學(xué)超聲基礎(chǔ)(北京:科學(xué)出版社)第415—418頁(yè)]

[29]Pennes H H 1948 J.Appl.Physiol.1 93

[30]Chen Z J,Zhou G L 2014 Inform.ComMun.4 36(in Chinese)[陳姝君,周廣麗 2014信息通信 4 36]

[31]Chen Z J 2008 Ph.D.D issertation(Nan jing:Nan jing University of Science and Technology)(in Chinese)[陳姝君2008博士學(xué)位論文(南京:南京理工大學(xué))]

[32]Li G Y 2011 M.S.D issertation(Beijing:Beijing Jiaotong University)(in Chinese)[黎光宇2011碩士學(xué)位論文(北京:北京交通大學(xué))]

[33]Xu G Z,Li Y,Yang S,Wu H L,Zhang S,Zhang J J 2010 BioMedical E lectrical IMpedance ToMograph y(Beijing:Machinery Industry Press)pp33–34(in Chinese)[徐桂芝,李穎,楊碩,吳煥麗,張帥,張劍軍2010生物醫(yī)學(xué)電阻抗成像技術(shù)(北京:機(jī)械工業(yè)出版社)第33—34頁(yè)]

[34]Lin MF 2014M.S.Dissertation(Nan jing:Nanjing University of Posts and TelecomMunications)(in Chinese)[林明鋒2014碩士學(xué)位論文(南京:南京郵電大學(xué))]

[35]Zhang L 2014 Electr.Design Eng.22 184(in Chinese)[張麗2014電子設(shè)計(jì)工程22 184]

[36]Q in S L 2000 CoMputat.Phys.17 314(in Chinese)[秦世倫2000計(jì)算物理17 314]

[37]Ma H,Wang G 2009 COMSOL Mu ltiphysics Basic Operation Guide and FAQ(Beijing:China ComMunications Press)(in Chinese)[馬慧,王剛 2009 COMSOL Mu ltiphysics基本操作指南和常見(jiàn)問(wèn)題解答 (北京:人民交通出版社)]

[38]Hu B,Jiang L X,Huang Y 2006 Tech.Acoust.25 613(in English)[胡兵,姜立新,黃瑛2006聲學(xué)技術(shù) 25 613]

[39]Li G 2013 M.S.D issertation(T ian jing:T ian jing University of Science and Technology)(in Chinese)[黎鴿2013碩士學(xué)位論文(天津:天津科技大學(xué))]

[40]Bessonova O,Wilkens V 2013 J.Acoust.Soc.AMer.134 4213

[41]Sun JM,Yu J,Guo X S,Zhang D 2013 Acta Phys.Sin.62 054301(in Chinese)[孫健明,于潔,郭霞生,章東2013物理學(xué)報(bào)62 054301]

[42]Chen T,Fan T,Zhang W,Q iu Y Y,Tu J,Guo X S,Zhang D 2014 J.Appl.Phys.115 114902

[43]Chen T 2014 Ph.D.D issertation(Nan jing:Nan jing University)(in Chinese)[陳濤 2014博士學(xué)位論文 (南京:南京大學(xué))]

PACS:43.80.Ev,43.35.YbDOI:10.7498/aps.66.164301

*Pro ject supported by the National Natural Science Foundation of China(G rant Nos.11474166,11604156),the Natu ral Science Foundation of Jiangsu Province,China(G rant No.BK 20161013),the Postdoctoral Science Foundation of China(G rant No.2016M591874)and the Priority AcadeMic PrograMDevelopMent of Jiangsu H igher Education Institu tions,China.

?Corresponding author.E-Mail:Maqingyu@njnu.edu.cn

N on invasive teMperatu re Mon itoring for high intensity focused u ltrasound therapy based on electrical iMpedance toMography?

Guo Ge-Pu Su Hui-Dan Ding He-Ping Ma Qing-Yu?
(School of Physics and Technology,Nanjing NorMal University,Nanjing 210023,China)

8 Ap ril 2017;revised Manuscrip t

7 June 2017)

As a neWtreatment modality With little thermal damage and feWcellmetastases to surrounding normal tissues,high intensity focused ultrasound(HIFU)therapy is considered to be one of theMost proMising technologies for tuMor therapy in the 21stcentury.However,noninvasive teMperatureMonitoring for the focal region exhibits great significance of p recise thermal dosage control in HIFU treatment.By combining electrical iMpedancemeasurement and HIFU,an electrical iMpedance toMography(EIT)based teMperature Monitoring Method using surface voltages is proposed to reconstruct the distribution of electrical conductivity inside the focal p lane on the basis of the teMperature dependent electrical iMpedance of tissues.In theoretical study,a coMp rehensive systeMof EIT measurement during HIFU therapy is established.With the consideration of acoustic absorp tion in viscous tissues,three-diMensional HelMholtz equation for HIFU is siMp lified into two-diMensional axisymMetric cylindrical coordinates,and the characteristics of teMperature rising in the focal region are derived using Pennes bio-heat transfer equation.Then,by introducing the teMperatureconductivity relation into tissues,the processing Methods for electrical field and surface voltage in the focal region are constructed With constant current injection froMtwo symmetrical electrodes.In simu lation study,by app lying the experimental paraMeters of the focused transducer,the distributions of acoustic p ressure and teMperature are simu lated at a fixed acoustic power,and then the corresponding distributions of conductivity in the focal p lane are achieved at diff erent treatment times for centric and eccentric focusing.Furthermore,With the simu lations of current density and electricalpotentialgenerated by the rotating current injection froM16 pairsof symMetricalelectrodes,32×32 voltages are detected by the 32 surface electrodes p laced around the focal p lane of theModel.In conductivity iMage reconstruction,themodified NeWton-Raphson(MNR)algorithMis eMp loyed to conduct iterative calcu lation.It shows that With the increase of HIFU treatMent tiMe,the electrical conductivity in the focal region increases accordingly and reaches a MaximuMvalue in the center due to the highest acoustic p ressure and theMost energy accumulation.It is p roved that not only the position of the focal center,but also the conductivity distribution inside the focal region can be restored accurately by the p roposed EIT based reconstruction algorithm.The favorable results deMonstrate the feasibility of teMperaturemonitoring during HIFU therapy,and also provide a neWmethod of evaluating the noninvasive effi cacy and controlling the dose based on electrical iMpedanceMeasureMents.

high intensity focused ultrasound,electrical impedance tomography,temperaturemonitoring,surface electrode voltage

10.7498/aps.66.164301

?國(guó)家自然科學(xué)基金(批準(zhǔn)號(hào):11474166,11604156)、江蘇省自然科學(xué)基金(批準(zhǔn)號(hào):BK 20161013)、國(guó)家博士后基金(批準(zhǔn)號(hào):2016M591874)和江蘇高校優(yōu)勢(shì)學(xué)科資助的課題.

?通信作者.E-Mail:Maqingyu@n jnu.edu.cn

?2017中國(guó)物理學(xué)會(huì)C h inese P hysica l Society

http://Wu lixb.iphy.ac.cn