時間反演聚焦經(jīng)顱磁聲電刺激仿真與實驗研究*

張 帥 黃明輝 焦立鵬 郭 梁 周振宇 徐桂芝

(1 河北工業(yè)大學(xué) 省部共建電工裝備可靠性與智能化國家重點實驗室 天津 300130)

(2 河北工業(yè)大學(xué) 天津市生物電工與智能健康重點實驗室 天津 300130)

0 引言

神經(jīng)調(diào)控技術(shù)是開展神經(jīng)科學(xué)研究與神經(jīng)系統(tǒng)疾病臨床治療的基本手段?；陔姶鸥袘?yīng)原理，通過采用電場或磁場等非侵入式方式刺激腦神經(jīng)組織，來達到調(diào)節(jié)腦功能的目的，包括經(jīng)顱磁刺激(Transcranial magnetic stimulation，TMS)[1]、經(jīng)顱直流電刺激(Transcranial direct current stimulation，tDCS)[2]和經(jīng)顱超聲刺激(Transcranial ultrasound stimulation，TUS)等方法[3]。經(jīng)顱磁聲電刺激(Transcranial magnetic acoustic electrical stimulation，TMAES)是一種新型無創(chuàng)腦神經(jīng)調(diào)控技術(shù)，將超聲波和靜磁場作用于神經(jīng)組織產(chǎn)生電流，從而對神經(jīng)組織進行刺激，刺激區(qū)域由超聲波在神經(jīng)組織的聚焦區(qū)域決定，同時靜磁場在組織中的衰減遠小于交變磁場，因此TMAES 是一種高空間分辨率的深部有效神經(jīng)刺激。

Norton[4]通過計算麥克斯韋方程組，理論上證明靜磁場和超聲場可以產(chǎn)生足夠強的感應(yīng)電流對大腦進行刺激。楊少華[5]對靜磁場中帶電凝膠受超聲作用而產(chǎn)生的電信號以及超聲回波信號做了頻譜分析，指出二者中心頻率一致，初步驗證了該方法的可行性。李慧雨等[6]使用超聲換能器向靜置于磁場中的導(dǎo)電樣本打入超聲，并測量沿不同方向的聲壓信號以及電信號，分析了施加聲場和感應(yīng)電場之間的一致性。Yuan 等[7]和劉世坤等[8]分別通過分析刺激后大鼠局部場電位信號的變化以及動物行為學(xué)實驗，證明了磁聲耦合電刺激方法可形成有效刺激。張帥等[9-10]通過建立包括刺激深度、刺激強度、聚焦面積、焦域長度以及電場強度梯度系數(shù)在內(nèi)的評價指標，對經(jīng)顱磁聲電刺激感應(yīng)電場聚焦性能進行了有效評價；并且基于H-H 神經(jīng)元模型，仿真了TMAES下不同超聲參數(shù)對神經(jīng)元放電模式的影響，對于具體神經(jīng)調(diào)控和得到最優(yōu)TMAES 參數(shù)提供理論指導(dǎo)。

相控陣超聲換能器可實現(xiàn)高空間分辨率、深層次的有效刺激，但是顱骨對聲場的散射和衍射嚴重影響了焦點區(qū)域的聲場[11]。骨骼的聲阻抗約是軟組織的10 倍[12]。Yuldashev 等[13]在陣列換能器聚焦模擬中發(fā)現(xiàn)了隨聲壓變化的焦點偏移現(xiàn)象。為了實現(xiàn)更加精準的超聲聚焦，F(xiàn)ink 等[14]將時間反演法引入了聲學(xué)領(lǐng)域，系統(tǒng)理論分析的同時進行了大量的水聲實驗，實現(xiàn)了在空間和時間上的同步聚焦。高翔等[15]對時間反演法進行了深入的研究并進行改進，應(yīng)用于流體、固體和分層介質(zhì)，進行裂紋檢測和目標定位。該方法在生物醫(yī)學(xué)和磁聲成像等方面也有廣泛應(yīng)用[16-17]。蘇暢等將時間反轉(zhuǎn)法用于經(jīng)顱超聲聚焦，證明時間反轉(zhuǎn)法相比于傳統(tǒng)的時間延時方法可以使聚焦更準確[18-20]。

現(xiàn)有經(jīng)顱磁聲電刺激相關(guān)文獻大多使用超聲換能器加聲準直器或透鏡約束方式進行神經(jīng)刺激，所建立的聲場和感應(yīng)電場的聚焦性仍有待提高。本文采用相控陣超聲換能器，通過在顱腦內(nèi)設(shè)置虛擬點聲源仿真模擬超聲傳播過程，記錄各陣元接收到超聲信號的時間，按照后到先發(fā)原則進行反演聚焦。相比于傳統(tǒng)相控聚焦，可以準確地在預(yù)設(shè)點實現(xiàn)聚焦。通過搭建實驗平臺，利用水聽器檢測到的數(shù)值繪制聚焦點處的歸一化聲壓圖，證實了時間反演法可以解決焦點偏移，實現(xiàn)自適應(yīng)聚焦，并對聲場和感應(yīng)電場分布一致性進行驗證，進一步證實了時間反演用于經(jīng)顱磁聲電刺激的可行性。

1 原理

1.1 TMAES原理

神經(jīng)組織樣本內(nèi)部帶電粒子受超聲波作用而振動，振動的粒子在磁場B0中切割磁感線產(chǎn)生與超聲方向垂直的洛倫茲力，導(dǎo)致正負電荷向相反方向移動，進而產(chǎn)生感應(yīng)電場E，經(jīng)顱磁聲電刺激的原理如圖1所示。超聲作用下樣本質(zhì)點沿y方向的振動速度為v(y，t)，q為樣本質(zhì)點所帶電荷量，電荷受到的洛倫茲力為[4]

圖1 TMAES 原理圖Fig.1 TMAES schematic

神經(jīng)組織樣本內(nèi)部產(chǎn)生的感應(yīng)電場強度E為

組織樣本質(zhì)點的振動速度與質(zhì)點處有效聲壓p(y)之間的關(guān)系為

式(3)中，ρ和c分別為神經(jīng)組織樣本密度、超聲波在神經(jīng)組織樣本的傳播速度。將式(3)帶入式(2)中，可得振動質(zhì)點處有效感應(yīng)電場強度為

1.2 傳統(tǒng)相控聚焦

經(jīng)顱磁聲電刺激需要聚焦聲束來產(chǎn)生感應(yīng)電場以進行神經(jīng)刺激，傳統(tǒng)相控聚焦方法是計算各陣元時間延時來使聲束偏轉(zhuǎn)聚焦，延時時間的計算取決于中心頻率、尺寸大小、焦點位置以及焦距，選取探頭中心為參考點，當陣元數(shù)N為奇數(shù)時，各陣元相對參考點的時間延時為

其中，c是超聲傳播速度，l為焦距，d為探頭陣元間距，θ為偏轉(zhuǎn)角度，n即第n個陣元(n=0，±1，···，±(N -1)/1)，當陣元數(shù)為偶數(shù)時，公式(5)中n則改為(n+0.5)。

1.3 時間反演法

為解決焦點偏移現(xiàn)象，本文采用時間反演法進行聲束聚焦，在預(yù)設(shè)焦點處設(shè)置虛擬點源，各陣元接收聲波，按照后到先發(fā)的原則，使得發(fā)射時各陣元擁有準確的延遲發(fā)射時間，可以精準地聚焦在理想位置，原理如圖2所示。利用互相關(guān)算法[21]可以準確地計算出各個陣元的時間延時，具體的實施步驟為：第n個陣元接收到的信號為Xn(t)，經(jīng)過時間反轉(zhuǎn)后，其對應(yīng)的信號為Xn(T -t)。將此陣元的信號重新定義Xnr(t)，在每一個反轉(zhuǎn)信號添加一個時間常量t0，各陣元時間反轉(zhuǎn)信號為

圖2 時間反演聚焦原理圖Fig.2 Schematic diagram of time reversal focusing

其中，Kn為各陣元的補償系數(shù)，各陣元激勵信號強度不一樣，一般陣元中心處權(quán)重大于邊緣處陣元的權(quán)重，可提高主瓣能量，減少旁瓣數(shù)量和抑制柵瓣。

2 數(shù)值模擬過程與結(jié)果

2.1 模型建立

本文仿真所用模型數(shù)據(jù)皆從二維顱腦圖像中獲得。其中，顱骨的厚度為5 mm 到15 mm 不等，大腦側(cè)面長152 mm，正面長131 mm，顱骨的聲學(xué)參數(shù)密度ρ1和聲速c1分別為1912 kg/m3、2300 m/s，大腦的聲學(xué)參數(shù)密度ρ2和聲速c2分別為1030 kg/m3、1550 m/s[22]。如圖3(a)所示，設(shè)置32 個線陣元模擬超聲換能器，與顱骨之間設(shè)置純水，各陣元中心頻率0.5 MHz，陣元寬度a為λ/2，陣元間距d為3λ/4，陣元長20 mm，虛擬點源置于換能器中心點正下方40 mm 處，四周設(shè)置軟聲場邊界模擬無限大介質(zhì)用來吸收出射波，防止反射的聲波影響計算精度。為了防止大腦神經(jīng)組織受連續(xù)脈沖波刺激產(chǎn)生熱積累現(xiàn)象引發(fā)不良反應(yīng)，所以發(fā)射周期性正弦波超聲脈沖群，每個脈沖群含有的基波周期數(shù)Nc為100，聲脈沖群重復(fù)頻率(Pulse repetition frequency，PRF)為1000 Hz，脈沖重復(fù)數(shù)(Pulse repetition number，PRB)為200[23]?？臻g網(wǎng)格步長dx= dy= dz= 0.4 mm，步進時間為0.5 μs，以滿足精度和穩(wěn)定性的要求。

2.2 點源發(fā)射以及顱內(nèi)磁場分布

利用圖3(a)仿真模型進行仿真，聲波向四周發(fā)射，圖3(b)為腦組織內(nèi)虛擬點聲源發(fā)射單正弦脈沖傳播過程中的二維截面圖。超聲經(jīng)大腦、顱骨和水域到達換能器各個陣元。由圖3(b)可見，顱骨對超聲的能量衰減、傳播路徑有較大的影響。由于顱骨聲阻抗較高，只有小部分的超聲能量穿過顱骨，并因折射使傳播路徑發(fā)生變化，部分超聲能量通過反射回到顱內(nèi)。

圖3 仿真模型與點源發(fā)射聲波Fig.3 The simulation model and point source emits sound waves

記錄超聲波到達換能器各陣元時的聲壓幅值以及時間并繪制歸一化聲壓圖，圖4為17號～32號陣元接收到的歸一化聲壓信號圖(將圖2中換能器陣列右邊第一個定義為1 號陣元)。圖4中左邊出現(xiàn)第一個峰值的曲線為17號陣元的歸一化聲壓圖，最后一條為32號陣元，因為在非均勻介質(zhì)中傳播路徑的不同，距離虛擬點源近的陣元相比于較遠的陣元接收到的聲壓幅值較高，聲波到達時間較短。記錄每個陣元的歸一化聲壓峰值所對應(yīng)的時間，并將32號陣元作為參考陣元，計算其余陣元與參考陣元的時間差，即為各個陣元的延遲發(fā)射時間。

圖4 部分陣元接收信號Fig.4 Some elements receive signal

永磁體提供經(jīng)顱磁聲電刺激所需的靜磁場，圖5為基于顱腦結(jié)構(gòu)的TMAES 系統(tǒng)模型，顱骨兩邊各放置一塊直徑80 mm、厚20 mm 的柱形永磁體。圖6為預(yù)設(shè)點源處沿永磁體中軸線上磁通密度分布曲線，曲線對稱分布，z軸0 mm 處為極小值0.5 T。將磁通密度降至最大值95%的區(qū)域近似看成磁場均勻分布區(qū)域，則分布區(qū)域為-12～12 mm，超聲聚焦區(qū)域的大小為毫米級別的，所以在聚焦點處的磁場可近似看成0.5 T 的均勻磁場，滿足磁聲電刺激的基本要求。

圖5 TMAES 系統(tǒng)模型Fig.5 TMAES system model

圖6 磁通密度模Fig.6 Norm of magnetic flux density

2.3 仿真結(jié)果分析

超聲換能器發(fā)射的聲波經(jīng)水域、顱骨和大腦在預(yù)設(shè)點上聚焦，圖7(a)為傳統(tǒng)相控聚焦聲場切面圖，圖中紅色部分為聚焦區(qū)域，藍色為非聚焦區(qū)域，單位為聲壓Pa。焦點形狀狹長不規(guī)則，選用焦點橫向和縱向聲壓峰值分別衰減到0.707倍時的區(qū)域為焦域，所得焦點縱向分辨率28.4 mm，橫向分辨率5.7 mm，實際焦點位置相比于預(yù)設(shè)焦點徑向偏移4 mm，軸向偏移2.5 mm。因為顱骨的厚度不均勻，使得在聚焦區(qū)域左邊有長約8 mm、寬約2 mm 的焦域，可對非聚焦區(qū)域形成刺激。圖7(b)為運用時間反演法的聚焦聲場示意圖，焦點形狀規(guī)則分布，呈梭形，焦點縱向分辨率23.2 mm，橫向分辨率4.6 mm，相比于傳統(tǒng)相控聚焦，聚焦區(qū)域減小，實際焦點位置相比于預(yù)設(shè)焦點徑向偏移1 mm，軸向偏移0.5 mm，實際焦點與預(yù)設(shè)焦點吻合度較高，說明時間反演法在分層的非均勻介質(zhì)中具有良好的聚焦效果。

圖7 傳統(tǒng)相控聚焦和時間反演聚焦Fig.7 The traditional phased array focusing and time reversal focusing

由公式(4)可得，當磁感應(yīng)強度恒定為0.5 T時，神經(jīng)組織密度和超聲傳播速度為定值，感應(yīng)電場強度數(shù)值變化規(guī)律與聲壓大小呈線性關(guān)系，聲場分布可近似為電場分布。圖8中，實線為運用時間反演法焦點徑向聲壓最大處軸向截線電場分布，虛線為傳統(tǒng)相控聚焦焦點徑向聲壓最大處軸向截線電場分布，兩種聚焦方式下測得的聲壓最大值分別為1.44 MPa 和1.32 MPa，經(jīng)公式(4)計算得，運用時間反演法刺激靶點處的感應(yīng)電場峰值為0.48 V/m，大于傳統(tǒng)相控聚焦的感應(yīng)電場峰值0.44 V/m，聚焦強度略有增加。傳統(tǒng)相控聚焦的聚焦區(qū)域與所設(shè)虛擬點源位置有一定的偏差，難以實現(xiàn)精準聚焦，而運用時間反演法的聚焦區(qū)域與虛擬點源的位置基本吻合，這是由于顱骨折射聲波使傳播路徑改變，進而導(dǎo)致焦點偏移，同時提高聲束聚焦精度和提高焦點處感應(yīng)電場峰值。

圖8 電場徑向分布曲線Fig.8 Radial electric field distribution curve

時間反演法同樣適用于多點聚焦，通過控制兩組或者多組陣元可以實現(xiàn)兩個或者多個焦點。多數(shù)神經(jīng)類疾病的產(chǎn)生并不是孤立的大腦區(qū)域異常的結(jié)果，因此研究不同區(qū)域的多焦點刺激同樣具有重要應(yīng)用價值[24]。圖9為運用時間反演生成的兩點聚焦聲場分布圖，將32 個陣元分成兩組，預(yù)設(shè)的兩個虛擬點源先后發(fā)射聲波，分別記錄兩組陣元對應(yīng)預(yù)設(shè)點源接收到的信號，各陣元接收到的信號幅值不同，通過調(diào)節(jié)各個陣元發(fā)射聲壓的幅值，可以改變各焦點處聲壓能量，因為顱骨厚度不均勻，在超聲傳播過程中對超聲能量的阻擋程度不同，傳播路徑也不同，導(dǎo)致兩個焦點的形狀和大小不一樣。

圖9 兩點聚焦Fig.9 Two-point focus

3 實驗驗證

3.1 實驗設(shè)計

本文設(shè)計仿體實驗來對比傳統(tǒng)相控聚焦和時間反演聚焦的效果，利用本實驗室Verasonics 開放式多通道超聲研究平臺實現(xiàn)焦點可視化，可以對刺激位置準確定位。如圖10所示，在水箱中放置由CT模型3D 打印而成的仿體顱骨，選擇與真實顱骨聲學(xué)特性相近的聚酯材料。在矢狀縫與冠狀縫交點正下方20 mm 處固定64 陣元相控陣超聲探頭(P4-2v探頭，Verasonics，美國)。在顱骨仿體內(nèi)部距離超聲探頭中心正上方40 mm 處放置水聽器(HNR-0500，ONDA，美國)，水聽器穩(wěn)定粘合于千分尺，用于檢測顱骨內(nèi)部的聲壓信號。因為水聽器的輸出信號十分微弱，所以需要輔助使用聲壓放大器(5660B PREAMP，OLYMPUS，日本)放大超聲聲壓信號，將60 dB 前置聲壓放大器連接在示波器(TDS3014，Tektronix，美國)上，探頭由Verasonics 系統(tǒng)控制各陣元時間延時發(fā)射連續(xù)超聲脈沖波，以1 mm 為單位在水平方向移動水聽器，并檢測聲壓。

圖10 實驗平臺Fig.10 Experiment platform

此外，又設(shè)計檢測實驗來考察聲場與感應(yīng)電場之間的關(guān)系。如圖11 所示，在靜磁場中放入與超聲換能器處于同一水平面的導(dǎo)電樣本，將金屬導(dǎo)線放置于導(dǎo)電樣本中，并與超聲換能器探頭表面相距40 mm，在相對位置不變的情況下，控制Verasonics系統(tǒng)超聲相控陣的焦距以2 mm 為單位從20 mm遞增至60 mm，金屬導(dǎo)線連接差分放大器放大測量到的電信號，使用萬用表記錄數(shù)值，并使用水聽器檢測與金屬導(dǎo)線處于同一區(qū)域的超聲信號，將測得的兩種信號歸一化處理并對比。

3.2 實驗結(jié)果

本實驗使用的超聲相控陣參數(shù)如下：陣元數(shù)N為64，陣元寬度a為0.25 mm，陣元自身長度h為14 mm，陣元間距d為0.3 mm，相控陣總長l為18.9 mm，如圖12 所示。當陣元間距d遠小于超聲波長λ時，波束指向性較差[25]，為得到形狀規(guī)則的焦點，設(shè)置各陣元基波頻率為1 MHz。預(yù)先設(shè)定這些參數(shù)進行仿真，以獲得陣列延遲傳輸所需時間。圖13(a)為仿真模式下使用時間反演法得到的陣元中心處縱向截面聲場分布圖，焦域縱向分辨率約為21.3 mm，橫向分辨率約為1.5 mm。圖13(b)為聚焦實驗中計算機屏幕顯示的GUI 控制界面，其中，弧形區(qū)域為顱骨輪廓，圖中的一個亮斑即為水聽器尖端，移動水聽器，使示波器顯示的電壓值最大，即可視為聚焦區(qū)域的中心點，此外在探頭發(fā)射聲波聚焦的過程中不考慮剪切波的影響[26]。

圖12 探頭組成示意圖Fig.12 Schematic diagram of probe composition

圖13 單點聚焦和實驗界面Fig.13 Single point focus and experimental interface

圖14中實線為仿真模式下y軸等于0 mm處橫向歸一化聲壓圖，由圖14 可見聲壓最大值位于z軸0 mm 處，基本沒有發(fā)生焦點偏移現(xiàn)象，兩種聚焦方式出現(xiàn)聲壓峰值的位置不同，反演聚焦時在z= 0 mm 處達到峰值，相控聚焦則在z=-1 mm處達到峰值，而用時間反演法根據(jù)各采樣點繪制的歸一化聲壓圖與仿真模式下繪制的聲壓圖分布大致相似，相關(guān)系數(shù)為0.823，證明了時間反演法能夠解決焦點偏移現(xiàn)象，實現(xiàn)精準聚焦。

圖14 兩種聚焦方式聲場對比Fig.14 Sound field comparison of two focusing methods

將實驗檢測到的聲場與電場數(shù)據(jù)歸一化處理并繪制如圖15所示的歸一化強度圖，兩條曲線分布基本重合，隨著焦距的增加，聲壓和電場的強度均有所增強，并在40 mm 處達到峰值，焦距超過40 mm后，強度均有所下降。通過計算得到兩條曲線的相關(guān)系數(shù)為0.913，據(jù)此可得，換能器聲場和磁聲耦合產(chǎn)生的感應(yīng)電場分布具有較高的一致性。

圖15 電場與聲場分布Fig.15 Electric field and sound field distribution

為了便于進行比較，將水聽器經(jīng)前置放大器連接到數(shù)字示波器顯示界面的波形(圖16)記錄下來，得到時域響應(yīng)圖，并將實驗數(shù)據(jù)導(dǎo)入到Matlab設(shè)置采樣頻率進行頻譜分析，圖17(a)為超聲信號頻譜分析圖，橫軸為頻率，縱軸表示能量大?。煌韺?dǎo)電凝膠兩端的電流連接示波器顯示波形并記錄，利用Matlab 進行頻譜分析，得到圖17(b)的電流信號的頻譜分析圖。電流波形頻譜分析要比超聲波形的頻譜分析效果要好，這是因為測量超聲時噪聲影響較大，使得采集到的信號頻帶較寬，但仍可得到大部分集中在1 MHz 左右，而測量電流信號時導(dǎo)電凝膠雜質(zhì)較少，使得信號較為純凈，相對頻帶窄得多，超聲基波頻率為1 MHz，所測得的電流信號是頻率為1 MHz 的交流信號。兩種波形頻譜分析得到的中心頻率一致的結(jié)論，可證明磁場和聲場耦合確實可以產(chǎn)生電流進行刺激。

圖16 測量超聲波形和測量電流波形Fig.16 Measuring ultrasonic wave shape and measuring current waveform

圖17 超聲波形頻譜分析和電流波形頻譜分析Fig.17 Ultrasonic shape spectrum analysis and current waveform spectrum analysis

此外，為進一步驗證檢測到的電信號確為聲場與磁場共同作用下產(chǎn)生，將導(dǎo)電樣本分別放置于靜磁場和聲場中，均未檢測到電信號的產(chǎn)生，此實驗聯(lián)合聲電場檢測實驗綜合證明了經(jīng)顱磁聲電刺激可行性。

4 結(jié)論

通過建立基于真實顱骨和大腦的物理場模型，分別運用傳統(tǒng)相控聚焦方法和時間反演法對大腦進行經(jīng)顱磁聲電刺激，并對結(jié)果進行分析。由于傳統(tǒng)相控陣聚焦方法只考慮聚焦點和各陣元的相對位置，并未對真實經(jīng)顱超聲的傳播過程進行計算，顱骨的形狀、厚度、輪廓引起的相位失真、能量衰減和折射效應(yīng)等未得到考量，運用時間反演法能夠在非均勻介質(zhì)中實現(xiàn)自適應(yīng)聚焦。仿真結(jié)果證明基于時間反演的相位調(diào)控能夠較好地解決因傳播介質(zhì)的非均勻性而導(dǎo)致的焦區(qū)偏移問題，提高了聚焦精度。

時間反演聚焦時各陣元的延遲信號可以通過計算機模擬來獲得，從而避免在人體組織植入聲源造成損傷。通過搭建實驗平臺，對兩種聚焦方式的聲場歸一化處理，驗證了時間反演法可以補償焦點偏移，并通過對磁聲耦合產(chǎn)生的電場與聲場強度歸一化處理，分析得出了它們之間存在較高的一致性。但本實驗僅考慮了顱骨對超聲的影響，后續(xù)需進一步考慮頭皮、白質(zhì)、灰質(zhì)對超聲聚焦的影響?；跁r間反演的經(jīng)顱磁聲電刺激擁有高空間分辨率以及有效精準刺激的優(yōu)點，有助推動精準神經(jīng)調(diào)控的發(fā)展。