本發(fā)明屬于超聲檢測(cè)與超聲成像技術(shù)領(lǐng)域，特別涉及一種適用于腦部聚焦超聲空化實(shí)時(shí)監(jiān)控的三維無(wú)源成像方法及系統(tǒng)。

背景技術(shù)：

空化效應(yīng)是指介質(zhì)中的空化核在一些外加物理場(chǎng)(例如超聲場(chǎng)、激光場(chǎng)、微波場(chǎng))的作用下產(chǎn)生的膨脹、收縮、塌陷的動(dòng)力學(xué)過(guò)程。腦部聚焦超聲治療主要是利用包膜微泡在超聲激勵(lì)下的空化效應(yīng)來(lái)達(dá)到治療腦部疾病的目的。已有研究表明，聚焦超聲聯(lián)合微泡開(kāi)放血腦屏障的主要機(jī)制是空化效應(yīng)，在較低聲壓下微泡發(fā)生穩(wěn)態(tài)空化，此時(shí)與大腦毛細(xì)血管尺寸相當(dāng)?shù)奈⑴菰谡駝?dòng)過(guò)程中與毛細(xì)血管相互作用，將構(gòu)成血腦屏障的內(nèi)皮細(xì)胞分開(kāi)，從而有利于血腦屏障的開(kāi)放。另外，空化效應(yīng)已被證實(shí)為超聲溶栓的主要機(jī)理，超聲空化效應(yīng)使栓塞部位的包膜微泡發(fā)生振動(dòng)及瞬間破裂，血栓表面變軟，通過(guò)損壞栓塊中的纖維蛋白結(jié)構(gòu)，有利于血栓的崩潰瓦解，增加溶栓藥物與栓塊的結(jié)合位點(diǎn)，從而加速溶栓。作為一種人工空化核，包膜微泡可以極大地降低空化閾值從而減少腦組織的額外損傷。然而，微泡在診斷超聲聲壓下的慣性收縮力也可能會(huì)導(dǎo)致微泡坍塌，從而對(duì)腦部細(xì)胞結(jié)構(gòu)產(chǎn)生破壞或者對(duì)血管壁形成不必要的損傷。由于空化效應(yīng)瞬時(shí)、隨機(jī)、不可控，因此對(duì)腦部聚焦超聲治療中的空化效應(yīng)進(jìn)行實(shí)時(shí)監(jiān)控是非常必要的。

近年來(lái)，國(guó)內(nèi)外研究表明，磁共振引導(dǎo)的經(jīng)顱聚焦超聲技術(shù)成為腦部超聲治療的新手段。磁共振監(jiān)控的優(yōu)勢(shì)主要體現(xiàn)在較高的圖像分辨率、對(duì)組織溫度實(shí)時(shí)變化的敏感性以及與聚焦超聲治療過(guò)程互不干擾等方面，然而常規(guī)的磁共振成像主要關(guān)注組織結(jié)構(gòu)的變化，難以實(shí)現(xiàn)對(duì)空化瞬態(tài)物理過(guò)程的實(shí)時(shí)檢測(cè)與監(jiān)控。傳統(tǒng)以聚焦波為發(fā)射方式的b模式超聲成像在腦部監(jiān)控存在超聲信號(hào)衰減、探測(cè)深度以及檢測(cè)靈敏度等相互制約的矛盾。且為了保證b模式成像與聚焦超聲信號(hào)互不影響，b模式成像只能在聚焦超聲脈沖發(fā)射之后或者不連續(xù)的聚焦超聲激勵(lì)下提供反饋圖像信息，無(wú)法對(duì)聚焦超聲作用過(guò)程中的空化進(jìn)行監(jiān)控，因此b模式成像也無(wú)法實(shí)現(xiàn)對(duì)腦部超聲治療過(guò)程中空化的實(shí)時(shí)監(jiān)控成像。以平面波為發(fā)射方式的超聲成像極大地提高了成像幀頻，但是也只能對(duì)聚焦超聲治療完成之后的空化進(jìn)行快速的監(jiān)控，依然無(wú)法與聚焦超聲治療信號(hào)做到同步。

因此，為了克服上述問(wèn)題，近年來(lái)興起了基于診斷超聲換能器的無(wú)源成像技術(shù)。無(wú)源成像技術(shù)是使診斷超聲換能器不發(fā)射脈沖信號(hào)只接收聚焦超聲換能器焦點(diǎn)處的空化信號(hào)，因此可以對(duì)聚焦超聲作用過(guò)程中的空化效應(yīng)進(jìn)行實(shí)時(shí)監(jiān)控且相比主動(dòng)超聲成像技術(shù)具有更高的檢測(cè)靈敏度。無(wú)源成像技術(shù)根據(jù)無(wú)源波束合成方法獲得空化源的空間分布，以此來(lái)對(duì)聚焦超聲治療進(jìn)行實(shí)時(shí)監(jiān)控。聚焦超聲治療監(jiān)控的準(zhǔn)確性取決于無(wú)源波束合成方法的性能。當(dāng)前，腦部聚焦超聲空化的無(wú)源成像技術(shù)主要是根據(jù)ct數(shù)據(jù)獲取成像區(qū)域有效聲速，然后使用基于延時(shí)-疊加-積分的時(shí)間曝光聲學(xué)(timeexposureacoustic,tea)方法進(jìn)行成像。由于受成像深度及換能器孔徑的限制，此方法的橫向/軸向分辨率較差，且微泡間的相互作用會(huì)形成嚴(yán)重的成像偽影。在非腦部聚焦超聲空化的無(wú)源成像中，魯棒的capon波束合成(robustcaponbeamformer,rcb)方法通過(guò)對(duì)每個(gè)成像位置的導(dǎo)向矢量重新估計(jì)并計(jì)算加權(quán)矢量，可以有效地消除微泡相互作用產(chǎn)生的偽影并同時(shí)提高橫向/軸向分辨率，但成像分辨率依然較差且成像偽影依然不能被完全消除。且當(dāng)前無(wú)源成像技術(shù)均基于一維線(xiàn)陣超聲換能器，只能提供腦部治療過(guò)程中空化效應(yīng)的平面分布，不能得到空化的三維空間分布，無(wú)法滿(mǎn)足腦部聚焦超聲治療實(shí)時(shí)監(jiān)控的臨床需求。

鑒于以上內(nèi)容，很有必要提出一種針對(duì)腦部聚焦超聲空化實(shí)時(shí)監(jiān)控的三維無(wú)源成像方法。

技術(shù)實(shí)現(xiàn)要素：

針對(duì)磁共振成像對(duì)腦部空化微泡信號(hào)的檢測(cè)不靈敏、傳統(tǒng)b模式和平面波主動(dòng)成像的探測(cè)深度和檢測(cè)靈敏度互相制約且無(wú)法監(jiān)控聚焦超聲作用過(guò)程中的空化效應(yīng)、以及現(xiàn)有二維無(wú)源成像技術(shù)在分辨率低、偽影高和無(wú)法提供空化三維空間信息等方面的問(wèn)題，本發(fā)明的目的在于提供一種用于腦部聚焦超聲空化實(shí)時(shí)監(jiān)控的三維無(wú)源成像方法及系統(tǒng)。

為實(shí)現(xiàn)上述目的，本發(fā)明采用了以下技術(shù)方案：

一種用于腦部聚焦超聲空化實(shí)時(shí)監(jiān)控的三維無(wú)源成像方法，包括以下步驟：

步驟一：計(jì)算超聲面陣換能器每個(gè)陣元方向上超聲波的有效聲速；

步驟二：利用拉格朗日插值法對(duì)超聲面陣換能器接收的聚焦超聲換能器焦點(diǎn)區(qū)域的空化信號(hào)進(jìn)行插值，然后根據(jù)步驟一所得有效聲速對(duì)空化信號(hào)進(jìn)行延時(shí)；

步驟三：構(gòu)造第j列陣元經(jīng)插值及延時(shí)后的空化信號(hào)的協(xié)方差矩陣，然后基于魯棒的capon波束合成方法計(jì)算最優(yōu)加權(quán)系數(shù)；對(duì)第j列陣元經(jīng)插值及延時(shí)后的空化信號(hào)進(jìn)行二值化并計(jì)算相位相干系數(shù)；根據(jù)該最優(yōu)加權(quán)系數(shù)以及相位相干系數(shù)，計(jì)算第j列陣元的空化源強(qiáng)度波束合成信號(hào)；

步驟四：根據(jù)步驟三計(jì)算的超聲面陣換能器各列陣元的空化源強(qiáng)度波束合成信號(hào)構(gòu)造協(xié)方差矩陣，然后基于魯棒的capon波束合成方法計(jì)算最優(yōu)加權(quán)系數(shù)；對(duì)各列陣元的空化源強(qiáng)度波束合成信號(hào)進(jìn)行二值化并計(jì)算相位相干系數(shù)；

步驟五：根據(jù)各列陣元的空化源強(qiáng)度波束合成信號(hào)以及步驟四所得最優(yōu)加權(quán)系數(shù)和相位相干系數(shù)計(jì)算超聲面陣換能器所有陣元的空化源強(qiáng)度波束合成信號(hào)和每個(gè)成像位置的空化源能量，得到空化源三維體數(shù)據(jù)，然后進(jìn)行腦部聚焦超聲空化的三維顯示。

上述步驟一中，每個(gè)陣元方向上超聲波的有效聲速按照以下公式計(jì)算：

其中k為除顱骨以外不同介質(zhì)的數(shù)目，li,j為空化信號(hào)到達(dá)超聲面陣換能器陣元(i,j)的距離，di,j,k和ci,j,k分別為空化信號(hào)到達(dá)陣元(i,j)的聲傳播路徑上除顱骨以外不同介質(zhì)的厚度和聲速，di,j,s為空化信號(hào)到達(dá)陣元(i,j)的聲傳播路徑上顱骨對(duì)應(yīng)位置的厚度，ci,j,s為超聲波在所述顱骨對(duì)應(yīng)位置的聲速；

超聲波在所述顱骨對(duì)應(yīng)位置的聲速ci,j,s根據(jù)以下公式計(jì)算：

ci,j,s＝c1φi,j,s+c2(1-φi,j,s)

其中φi,j,s為所述顱骨對(duì)應(yīng)位置的孔隙度，c1為超聲波在水中的傳播速度，c2為超聲波在顱骨中的最大傳播速度。

上述步驟二，具體包括以下步驟：

(2.1)設(shè)置超聲面陣換能器的陣元變跡，使換能器工作在不發(fā)射脈沖信號(hào)只接收信號(hào)的模式，同步觸發(fā)聚焦超聲換能器和超聲面陣換能器，由超聲面陣換能器采集聚焦超聲換能器作用過(guò)程中的焦點(diǎn)處空化的三維射頻信號(hào)pi,j，其中i＝1,2,...,m；j＝1,2,...,n；m和n分別為超聲面陣換能器y和x方向的陣元數(shù)目；

(2.2)根據(jù)以下拉格朗日插值多項(xiàng)式對(duì)每個(gè)陣元采集到的射頻信號(hào)進(jìn)行插值：

其中，為插值后的空化信號(hào)，ki,j(tl,t)為拉格朗日基本多項(xiàng)式，pi,j(tl)為空化信號(hào)采樣點(diǎn)，p為采樣點(diǎn)個(gè)數(shù)；

(2.3)根據(jù)三維成像區(qū)域和成像精度，對(duì)步驟(2.2)中所得信號(hào)進(jìn)行延時(shí)處理，得到延時(shí)后的信號(hào)si,j(x,y,z,t)：

其中di,j(x,y,z)為成像位置(x,y,z)到超聲面陣換能器陣元位置的空間距離，ci,j為陣元方向上超聲波的有效聲速。

上述步驟三，具體包括以下步驟：

(3.1)在一定的時(shí)間區(qū)間[t0,t0+δt]內(nèi)構(gòu)造協(xié)方差矩陣：

其中sj(x,y,z,t)為第j列陣元經(jīng)插值及延時(shí)后的空化信號(hào)矩陣，[·]^t代表矩陣的轉(zhuǎn)置；

(3.2)對(duì)步驟(3.1)所得協(xié)方差矩陣做特征值分解：

其中vj為特征值矩陣，uj中第i列為vj中對(duì)角線(xiàn)上第i個(gè)元素對(duì)應(yīng)的特征向量；

(3.3)根據(jù)步驟(3.2)中所得特征值矩陣vj和特征向量矩陣uj以及拉格朗日迭代法計(jì)算每個(gè)成像位置的導(dǎo)向矢量

其中為m行1列(m×1)預(yù)設(shè)導(dǎo)向矢量，λ為拉格朗日乘數(shù)，i為m行m列(m×m)單位矩陣；

(3.4)基于魯棒的capon波束合成方法并利用步驟(3.1)所得協(xié)方差矩陣rj和步驟(3.3)所得導(dǎo)向矢量計(jì)算最優(yōu)加權(quán)系數(shù)wj：

其中||·||代表歐幾里得范數(shù)，[·]^-1代表矩陣的逆；

(3.5)將第j列陣元經(jīng)插值及延時(shí)后的空化信號(hào)si,j(x,y,z,t)中大于零的點(diǎn)置為1，小于等于零的點(diǎn)置為-1，得到二值化信號(hào)bi,j(x,y,z,t)，然后計(jì)算相位相干系數(shù)

其中p為可調(diào)參數(shù)；

(3.6)利用步驟(3.4)所得最優(yōu)加權(quán)系數(shù)以及步驟(3.5)所得相位相干系數(shù)計(jì)算第j列陣元的空化源強(qiáng)度波束合成信號(hào)：

上述步驟四，具體包括以下步驟：

(4.1)計(jì)算超聲面陣換能器各列陣元的空化源強(qiáng)度波束合成信號(hào)矩陣：

q(x,y,z,t)＝[q1(x,y,z,t)；q2(x,y,z,t)；...；qn(x,y,z,t)]

其中qj(x,y,z,t)為第j列陣元的空化源強(qiáng)度波束合成信號(hào)；

(4.2)在一定的時(shí)間區(qū)間[t0,t0+δt]內(nèi)，根據(jù)步驟(4.1)所得空化源強(qiáng)度波束合成信號(hào)矩陣q(x,y,z,t)計(jì)算協(xié)方差矩陣：

其中[·]^t代表矩陣的轉(zhuǎn)置；

(4.3)對(duì)步驟(4.2)所得協(xié)方差矩陣做特征值分解：

r(x,y,z)＝uvu^t

其中v為特征值矩陣，u中第j列為v中對(duì)角線(xiàn)上第j個(gè)元素對(duì)應(yīng)的特征向量，j＝1,2,...,n；

(4.4)根據(jù)步驟(4.3)中所得特征值矩陣v和特征向量矩陣u以及拉格朗日迭代法計(jì)算每個(gè)成像位置的導(dǎo)向矢量

其中為n行1列(n×1)預(yù)設(shè)導(dǎo)向矢量，λ為拉格朗日乘數(shù)，i為n行n列(n×n)單位矩陣；

(4.5)基于魯棒的capon波束合成方法并利用步驟(4.2)所得協(xié)方差矩陣r和步驟(4.4)所得導(dǎo)向矢量計(jì)算最優(yōu)加權(quán)系數(shù)w：

其中||·||代表歐幾里得范數(shù)，[·]^-1代表矩陣的逆；

(4.6)將每列陣元的空化源強(qiáng)度波束合成信號(hào)qj(x,y,z,t)中大于零的點(diǎn)置為1，小于等于零的點(diǎn)置為-1，得到二值化信號(hào)bj(x,y,z,t)，然后計(jì)算相位相干系數(shù)scf^p(x,y,z,t)：

其中p為可調(diào)參數(shù)。

上述步驟五中，所有陣元的空化源強(qiáng)度波束合成信號(hào)按照以下公式計(jì)算：

q(x,y,z,t)＝w^tq(x,y,z,t)scf^p(x,y,z,t)

在一定的時(shí)間區(qū)間[t0,t0+δt]內(nèi)對(duì)q(x,y,z,t)的平方進(jìn)行積分，得到每個(gè)成像位置處的空化源能量i(x,y,z)：

上述步驟五中，三維顯示具體步驟為：

(5.1)對(duì)所得的三維體數(shù)據(jù)進(jìn)行閾值化處理，閾值一般取體數(shù)據(jù)最大值的0.1～0.2倍；

(5.2)選擇濾波器對(duì)步驟(5.1)所得閾值化后的三維體數(shù)據(jù)進(jìn)行平滑處理，然后提取出體素等值面；

(5.3)根據(jù)步驟(5.2)所得等值面構(gòu)造三維曲面，然后設(shè)置顏色、光照、視角等進(jìn)行三維顯示。

一種用于腦部聚焦超聲空化實(shí)時(shí)監(jiān)控的三維無(wú)源成像系統(tǒng)，包括超聲面陣換能器、聲速計(jì)算模塊、空化信號(hào)校準(zhǔn)模塊、自適應(yīng)波束合成模塊以及三維顯示模塊；

所述聲速計(jì)算模塊計(jì)算超聲面陣換能器每個(gè)陣元方向上超聲波的有效聲速(參見(jiàn)上述步驟一)；

所述空化信號(hào)校準(zhǔn)模塊對(duì)超聲面陣換能器接收的聚焦超聲換能器焦點(diǎn)區(qū)域的空化信號(hào)進(jìn)行插值以及延時(shí)處理(參見(jiàn)上述步驟二)；

所述自適應(yīng)波束合成模塊對(duì)超聲面陣換能器兩個(gè)方向的陣元做魯棒的capon波束合成，并基于空化信號(hào)相位間的差異引入相位相干系數(shù)對(duì)波束合成算法進(jìn)行修正；

所述三維顯示模塊利用自適應(yīng)波束合成模塊得到的空化源三維體數(shù)據(jù)進(jìn)行腦部聚焦超聲空化的三維顯示。

本發(fā)明的有益效果體現(xiàn)在：

本發(fā)明是一種針對(duì)腦部聚焦超聲空化實(shí)時(shí)監(jiān)控的三維無(wú)源成像技術(shù)，針對(duì)磁共振成像、傳統(tǒng)b模式和平面波主動(dòng)成像以及現(xiàn)有無(wú)源成像方法在監(jiān)控腦部聚焦超聲治療方面存在的缺陷，通過(guò)同步觸發(fā)聚焦超聲治療換能器和超聲面陣換能器采集治療過(guò)程中空化的三維射頻信號(hào)并根據(jù)ct掃描數(shù)據(jù)對(duì)信號(hào)延時(shí)進(jìn)行校準(zhǔn)?；隰敯舻腸apon波束合成方法對(duì)面陣陣元的延時(shí)信號(hào)進(jìn)行波束合成，對(duì)來(lái)自其他方向的信號(hào)進(jìn)行了大幅度的抑制；根據(jù)通道信號(hào)相位間的差異引入了相位相干系數(shù)對(duì)波束合成算法進(jìn)行改進(jìn)，提高了成像分辨率，從而準(zhǔn)確地反映了空化在三維空間中的分布，對(duì)于腦部聚焦超聲治療過(guò)程中空化效應(yīng)的實(shí)時(shí)監(jiān)控有著重要的意義。

附圖說(shuō)明

圖1為腦部ct掃描數(shù)據(jù)及用于三維無(wú)源成像方法的顱骨內(nèi)聲速獲取示意圖；

圖2為超聲面陣換能器獲取的空化信號(hào)的插值及延時(shí)流程圖；

圖3為某超聲面陣換能器有無(wú)顱骨遮擋時(shí)空化信號(hào)的時(shí)頻域?qū)Ρ葓D；

圖4為超聲面陣換能器第j列陣元的空化源強(qiáng)度波束合成的流程圖；

圖5為超聲面陣換能器每列陣元波束合成中所得的協(xié)方差矩陣、特征向量矩陣、特征值、導(dǎo)向矢量、最優(yōu)加權(quán)系數(shù)和相位相干系數(shù)結(jié)果；

圖6為超聲面陣換能器所有陣元的空化源強(qiáng)度波束合成的流程圖；

圖7為超聲面陣換能器所有陣元波束合成中所得的協(xié)方差矩陣、特征向量矩陣、特征值、導(dǎo)向矢量、最優(yōu)加權(quán)系數(shù)和相位相干系數(shù)結(jié)果；

圖8為基于不同無(wú)源成像算法的三維空化成像x-z切面結(jié)果；

圖9為三維無(wú)源成像顯示流程圖；

圖10為不同條件下的三維無(wú)源成像顯示結(jié)果圖。

具體實(shí)施方式

下面結(jié)合附圖和實(shí)施例對(duì)本發(fā)明做詳細(xì)描述。

本發(fā)明涉及到的用于腦部聚焦超聲空化實(shí)時(shí)監(jiān)控的三維無(wú)源成像方法包括以下具體步驟：

步驟一：根據(jù)顱骨ct掃描數(shù)據(jù)和基于顱骨孔隙度的經(jīng)驗(yàn)?zāi)Ｐ陀?jì)算超聲波在顱骨中的傳播速度，計(jì)算超聲面陣換能器每個(gè)陣元方向上超聲波的有效聲速。其具體步驟如下：

(1.1)獲取顱骨ct掃描切片數(shù)據(jù)，圖1中(a)為獲取的顱骨ct掃描切片數(shù)據(jù)，(b)為利用(a)數(shù)據(jù)得到的顱骨三維重建結(jié)果。根據(jù)圖1(a)所示的顱骨ct掃描切片數(shù)據(jù)計(jì)算空化信號(hào)到達(dá)陣元(i,j)的聲傳播路徑上顱骨對(duì)應(yīng)位置的孔隙度：

其中hi,j為亨氏單位；

(1.2)在顱骨孔隙度經(jīng)驗(yàn)?zāi)Ｐ偷幕A(chǔ)上計(jì)算超聲波在顱骨對(duì)應(yīng)位置的聲速ci,j,s：

ci,j,s＝c1φi,j,s+c2(1-φi,j,s)

其中c1為超聲波在水中的傳播速度，c2為超聲波在顱骨中的最大傳播速度；

(1.3)圖1中(c)為顱骨孔隙度計(jì)算結(jié)果，(d)和(e)分別為顱骨內(nèi)部密度和超聲波傳播速度的分布結(jié)果，密度和聲速分布圖像基本一致，證實(shí)超聲波在密度較大的地方傳播較快；且兩幅圖像均顯示顱骨內(nèi)側(cè)和外層密度和聲速較大，而中間層減?。贿@是因?yàn)轱B骨內(nèi)外兩側(cè)為密質(zhì)層而中間層為松質(zhì)層，有靜脈通過(guò)。(f)為(e)的局部放大結(jié)果，用于三維無(wú)源成像中有效聲速的計(jì)算；成像區(qū)域選擇顱骨較薄的顳窗部位，以盡可能減小空化信號(hào)的衰減；

(1.4)在實(shí)際實(shí)驗(yàn)過(guò)程中，通過(guò)測(cè)量空化信號(hào)到超聲面陣換能器傳播路徑上不同介質(zhì)(例如顱骨、液體、腦組織、腦組織仿體等)的厚度和聲速，計(jì)算每個(gè)陣元方向上的有效聲速

其中k為除顱骨以外不同介質(zhì)的數(shù)目，li,j為空化信號(hào)到達(dá)超聲面陣換能器陣元(i,j)的距離，di,j,k和ci,j,k分別為空化信號(hào)到達(dá)陣元(i,j)聲傳播路徑上除顱骨以外不同介質(zhì)的厚度和聲速，di,j,s為顱骨對(duì)應(yīng)位置的厚度，ci,j,s為根據(jù)步驟(1.2)計(jì)算的聲速。

步驟二：超聲面陣換能器接收聚焦超聲換能器焦點(diǎn)區(qū)域的空化信號(hào)，利用拉格朗日插值法對(duì)空化信號(hào)進(jìn)行插值，并根據(jù)步驟一所得有效聲速對(duì)空化信號(hào)進(jìn)行延時(shí)，得到延時(shí)信號(hào)。其具體步驟如下(圖2)：

(2.1)設(shè)置超聲面陣換能器的陣元變跡，使換能器工作在不發(fā)射脈沖信號(hào)只接收信號(hào)的模式，通過(guò)波形發(fā)生器編寫(xiě)觸發(fā)波形來(lái)同步觸發(fā)聚焦超聲換能器和面陣換能器，采集聚焦超聲換能器作用過(guò)程中的焦點(diǎn)處空化的三維射頻信號(hào)pi,j，其中i＝1,2,...,m；j＝1,2,...,n；m和n分別為面陣換能器y和x方向的陣元數(shù)目(即面陣中陣元的行數(shù)和列數(shù))；

(2.2)由于實(shí)際實(shí)驗(yàn)中信號(hào)采樣率的限制，可能會(huì)造成信號(hào)延時(shí)的精度不夠，因此在信號(hào)延時(shí)之前對(duì)每個(gè)陣元采集到的空化信號(hào)進(jìn)行插值；假設(shè)原信號(hào)采樣點(diǎn)為pi,j(t0),pi,j(t1),...,pi,j(tp-1)，對(duì)p個(gè)采樣點(diǎn)進(jìn)行拉格朗日插值，計(jì)算當(dāng)前列當(dāng)前陣元的拉格朗日基本多項(xiàng)式：

(2.3)設(shè)置插值精度，一般在每?jī)蓚€(gè)采樣點(diǎn)之間插入4～8個(gè)點(diǎn)，根據(jù)拉格朗日插值多項(xiàng)式計(jì)算原采樣點(diǎn)之外的插值點(diǎn)的函數(shù)值，當(dāng)前列當(dāng)前陣元的拉格朗日插值多項(xiàng)式為：

其中ki,j(tl,t)為拉格朗日基本多項(xiàng)式，為插值后的空化信號(hào)；

(2.4)跳轉(zhuǎn)至當(dāng)前列下一個(gè)陣元重復(fù)步驟(2.2)和(2.3)，直到當(dāng)前列所有陣元采集到的空化信號(hào)均插值完畢；

(2.5)跳轉(zhuǎn)至下一列陣元并重復(fù)步驟(2.2)～(2.4)，直到面陣所有陣元采集到的空化信號(hào)均插值完畢；

(2.6)規(guī)劃三維成像區(qū)域和成像精度，一般三維成像區(qū)域橫向選擇為面陣換能器的尺寸，軸向選擇為成像區(qū)域與面陣之間距離的2～3倍；當(dāng)需要觀(guān)察局部放大的空化分布時(shí)，可適當(dāng)?shù)販p小成像區(qū)域。成像精度一般為所選擇的成像區(qū)域中有(100～200)×(100～200)×(100～200)個(gè)體素，成像精度根據(jù)需要可以適當(dāng)增加或減小；

(2.7)計(jì)算每個(gè)成像位置(x,y,z)到面陣換能器當(dāng)前列當(dāng)前陣元位置(xi,yj,0)的空間距離：

(2.8)結(jié)合步驟(1.4)計(jì)算得到的有效聲速ci,j對(duì)插值后的空化信號(hào)做延時(shí)處理，延時(shí)后的信號(hào)為：

(2.9)跳轉(zhuǎn)至當(dāng)前列下一個(gè)陣元并重復(fù)步驟(2.7)和(2.8)，直到當(dāng)前列所有陣元采集到的信號(hào)均延時(shí)完畢；

(2.10)跳轉(zhuǎn)至下一列陣元并重復(fù)步驟(2.7)～(2.9)，直到面陣所有陣元采集到的信號(hào)均延時(shí)完畢；

如圖3所示，其中(a)為超聲面陣換能器陣列示意圖，由m×n個(gè)陣元組成，陣元均勻分布在32行32列柵格上，陣元總數(shù)為1024，孔徑大小為38.4×38.4mm。(b)和(c)為無(wú)顱骨遮擋條件下經(jīng)插值及延時(shí)的空化信號(hào)和相應(yīng)的頻譜分布；(d)和(e)為有顱骨遮擋條件下經(jīng)插值及延時(shí)的空化信號(hào)和相應(yīng)的頻譜分布。結(jié)果顯示有顱骨遮擋的空化信號(hào)強(qiáng)度明顯低于無(wú)顱骨遮擋情況，說(shuō)明顱骨對(duì)空化信號(hào)有一定的衰減；無(wú)顱骨遮擋下的空化信號(hào)頻譜成分尤其高頻成分豐富而有顱骨遮擋下空化信號(hào)中約7mhz以上的高頻成分被抑制，說(shuō)明顱骨對(duì)于空化信號(hào)相當(dāng)于一個(gè)低通濾波器。

步驟三：構(gòu)造每列陣元經(jīng)插值及延時(shí)后的空化信號(hào)的協(xié)方差矩陣，然后基于魯棒的capon波束合成方法(lij,stoicap,wangz.onrobustcaponbeamforminganddiagonalloading[j].signalprocessingieeetransactionson,2003,51(7):1702-1715.)計(jì)算最優(yōu)加權(quán)系數(shù)；對(duì)每列陣元經(jīng)插值及延時(shí)后的空化信號(hào)進(jìn)行二值化并計(jì)算相位相干系數(shù)；根據(jù)該最優(yōu)加權(quán)系數(shù)以及相位相干系數(shù)，計(jì)算每列陣元的空化源強(qiáng)度波束合成信號(hào)。其具體步驟如下(圖4)：

(3.1)提取第j列m個(gè)陣元經(jīng)過(guò)拉格朗日插值及延時(shí)之后的空化信號(hào)；

(3.2)在一定的時(shí)間區(qū)間[t0,t0+δt]內(nèi)構(gòu)造協(xié)方差矩陣：

其中sj(x,y,z,t)＝[s1,j(x,y,z,t)；s2,j(x,y,z,t)；...；sm,j(x,y,z,t)]為第j列m個(gè)陣元經(jīng)插值及延時(shí)后的空化信號(hào)矩陣，[·]^t代表矩陣的轉(zhuǎn)置；時(shí)間長(zhǎng)度δt取決于聚焦超聲每次照射的時(shí)間；

(3.3)對(duì)步驟(3.2)所得協(xié)方差矩陣做特征值分解：

其中vj為特征值矩陣，uj中第i列為vj中對(duì)角線(xiàn)上第i個(gè)元素(特征值)對(duì)應(yīng)的特征向量，i＝1,2,...,m；

(3.4)設(shè)置預(yù)設(shè)導(dǎo)向矢量和拉格朗日乘數(shù)λ，結(jié)合步驟(3.3)中所得特征值矩陣vj和特征向量矩陣uj計(jì)算每個(gè)成像位置的導(dǎo)向矢量

其中，i為單位矩陣(m×m)；

(3.5)基于魯棒的capon波束合成方法并利用步驟(3.2)所得協(xié)方差矩陣rj和步驟(3.4)所得導(dǎo)向矢量計(jì)算最優(yōu)加權(quán)系數(shù)wj：

其中||·||代表歐幾里得范數(shù)，[·]^-1代表矩陣的逆；

(3.6)將第j列當(dāng)前陣元經(jīng)插值及延時(shí)后的空化信號(hào)中大于零的點(diǎn)置為1，小于等于零的點(diǎn)置為-1，得到二值化信號(hào)：

(3.7)跳轉(zhuǎn)至下一個(gè)陣元并重復(fù)步驟(3.6)，直到第j列所有陣元經(jīng)插值及延時(shí)后的空化信號(hào)均已被二值化；

(3.8)根據(jù)步驟(3.7)所得第j列m個(gè)陣元的二值化信號(hào)計(jì)算相位相干系數(shù)

其中p為可調(diào)參數(shù)以改善圖像質(zhì)量，p的范圍一般為0～2；

(3.9)利用步驟(3.5)所得最優(yōu)加權(quán)系數(shù)以及步驟(3.8)所得相位相干系數(shù)計(jì)算第j列m個(gè)陣元的空化源強(qiáng)度波束合成信號(hào)：

如圖5所示，(a)為每列陣元經(jīng)插值及延時(shí)后的空化信號(hào)矩陣sj(x,y,z,t)所得到的協(xié)方差矩陣，(b)和(c)分別為(a)中協(xié)方差矩陣經(jīng)過(guò)特征值分解得到的特征向量矩陣和特征值，(d)為根據(jù)(b)中特征向量矩陣和(c)中特征值計(jì)算得到的導(dǎo)向矢量，(e)為利用(a)中協(xié)方差矩陣和(d)中導(dǎo)向矢量得到的最優(yōu)加權(quán)系數(shù)，(f)為利用步驟(3.8)得到的相位相干系數(shù)。

步驟四：根據(jù)步驟三計(jì)算的超聲面陣換能器各列陣元的空化源強(qiáng)度波束合成信號(hào)構(gòu)造協(xié)方差矩陣，然后基于魯棒的capon波束合成方法計(jì)算最優(yōu)加權(quán)系數(shù)；對(duì)各列陣元的空化源強(qiáng)度波束合成信號(hào)進(jìn)行二值化并計(jì)算相位相干系數(shù)。其具體步驟如下(圖6)：

(4.1)重復(fù)步驟(3.1)～(3.9)，直到第1～n列陣元的空化源強(qiáng)度波束合成信號(hào)均被計(jì)算完畢；

(4.2)構(gòu)造第1～n列陣元的空化源強(qiáng)度波束合成信號(hào)矩陣：

q(x,y,z,t)＝[q1(x,y,z,t)；q2(x,y,z,t)；...；qn(x,y,z,t)]

(4.3)根據(jù)步驟(4.2)所得第1～n列陣元的空化源強(qiáng)度波束合成信號(hào)矩陣q(x,y,z,t)計(jì)算協(xié)方差矩陣：

其中[·]^t代表矩陣的轉(zhuǎn)置；

(4.4)對(duì)步驟(4.3)所得協(xié)方差矩陣做特征值分解：

r(x,y,z)＝uvu^t

其中v為特征值矩陣，u中第j列為v中對(duì)角線(xiàn)上第j個(gè)元素(特征值)對(duì)應(yīng)的特征向量，j＝1,2,...,n；

(4.5)設(shè)置預(yù)設(shè)導(dǎo)向矢量和拉格朗日乘數(shù)λ，結(jié)合步驟(4.4)中所得特征值矩陣v和特征向量矩陣u計(jì)算每個(gè)成像位置的導(dǎo)向矢量

其中，i為單位矩陣(n×n)；

(4.6)基于魯棒的capon波束合成方法并利用步驟(4.3)所得協(xié)方差矩陣r和步驟(4.5)所得導(dǎo)向矢量計(jì)算最優(yōu)加權(quán)系數(shù)w：

其中||·||代表歐幾里得范數(shù)，[·]^-1代表矩陣的逆；

(4.7)將第j列陣元的空化源強(qiáng)度波束合成信號(hào)中大于零的點(diǎn)置為1，小于等于零的點(diǎn)置為-1，得到二值化信號(hào)bj(x,y,z,t)：

(4.8)跳轉(zhuǎn)至下一列陣元并重復(fù)步驟(4.7)，直到n列陣元的空化源強(qiáng)度波束合成信號(hào)均已被二值化；

(4.9)根據(jù)步驟(4.8)所得n個(gè)二值化信號(hào)計(jì)算相位相干系數(shù)scf^p(x,y,z,t)：

其中p為可調(diào)參數(shù)以改善圖像質(zhì)量，p的范圍一般為0～2；

如圖7所示，(a)為根據(jù)n列陣元的空化源強(qiáng)度波束合成信號(hào)矩陣q(x,y,z,t)所得到的協(xié)方差矩陣，(b)和(c)為(a)中協(xié)方差矩陣經(jīng)過(guò)特征值分解得到的特征向量矩陣和特征值，(d)為根據(jù)(b)中特征向量矩陣和(c)中特征值計(jì)算得到的導(dǎo)向矢量，(e)為利用(a)中協(xié)方差矩陣和(d)中導(dǎo)向矢量得到的最優(yōu)加權(quán)系數(shù)，(f)為利用步驟(4.9)得到的相位相干系數(shù)。

步驟五：根據(jù)步驟四所得每列陣元的空化源強(qiáng)度波束合成信號(hào)、最優(yōu)加權(quán)系數(shù)和相位相干系數(shù)計(jì)算所有陣元的空化源強(qiáng)度波束合成信號(hào)和每個(gè)成像位置的空化源能量，得到空化源三維體數(shù)據(jù)，對(duì)三維體數(shù)據(jù)進(jìn)行閾值化、平滑等處理，然后進(jìn)行腦部聚焦超聲空化的三維顯示。其具體步驟如下：

(5.1)利用步驟(4.2)所得第1～n列陣元的空化源強(qiáng)度波束合成信號(hào)矩陣q(x,y,z,t)、步驟(4.6)所得最優(yōu)加權(quán)系數(shù)以及步驟(4.9)所得相位相干系數(shù)計(jì)算所有陣元的空化源強(qiáng)度波束合成信號(hào)：

q(x,y,z,t)＝w^tq(x,y,z,t)scf^p(x,y,z,t)

(5.2)在一定的時(shí)間區(qū)間[t0,t0+δt]內(nèi)對(duì)步驟(5.1)所得q(x,y,z,t)的平方進(jìn)行積分，得到每個(gè)成像位置處的空化源能量(空化源三維體數(shù)據(jù))：

如圖8所示，(a)～(c)分別為使用傳統(tǒng)的tea和rcb算法以及本發(fā)明所述方法得到的三維無(wú)源成像的x-z切面，結(jié)果顯示傳統(tǒng)tea算法分辨率較差，成像偽影嚴(yán)重干擾了空化區(qū)域，rcb算法在一定程度上提高了分辨率并抑制了偽影，但在換能器遠(yuǎn)場(chǎng)仍存在較大的偽影，本發(fā)明則能在提高分辨率同時(shí)對(duì)遠(yuǎn)場(chǎng)的偽影進(jìn)行大幅度抑制。(d)和(e)分別為x、z方向上的分辨率曲線(xiàn)(虛線(xiàn)為tea，點(diǎn)虛線(xiàn)為rcb，實(shí)線(xiàn)為本發(fā)明方法)，(f)為偽影抑制程度的定量化(1～3分別代表tea，rcb和本發(fā)明方法)，定量結(jié)果顯示本發(fā)明所述方法具有良好的分辨性能，使得空化被局限在一個(gè)小區(qū)域內(nèi)，從而實(shí)現(xiàn)空化在三維空間的精準(zhǔn)定位。

如圖9所示，三維無(wú)源成像的顯示流程為：

(5.3)從步驟(5.2)獲取空化源三維體數(shù)據(jù)之后，設(shè)置閾值為所有體數(shù)據(jù)最大值的0.1～0.2倍，將體數(shù)據(jù)值小于此閾值的置為零；

(5.4)選擇濾波器為savitzky-golay濾波器，設(shè)置savitzky-golay濾波器的窗寬、多項(xiàng)式模型的階數(shù)以及平滑點(diǎn)數(shù)對(duì)步驟(5.3)所得閾值化后的體數(shù)據(jù)進(jìn)行平滑處理；

(5.5)計(jì)算步驟(5.4)所得平滑后的體數(shù)據(jù)的等值面的表面和頂點(diǎn)，提取出所有體素的等值面；

(5.6)根據(jù)步驟(5.5)所得等值面構(gòu)造三維曲面，計(jì)算等值面的頂點(diǎn)法線(xiàn)，最后設(shè)置顏色、光照、視角等進(jìn)行三維顯示。

如圖10所示，(a)～(d)分別為對(duì)不同空化微泡數(shù)目(20,50,100,200)的三維無(wú)源成像顯示結(jié)果，結(jié)果說(shuō)明本發(fā)明所述方法得到的三維無(wú)源成像具有足夠高的精度。

本發(fā)明以相位相干系數(shù)修正的魯棒的capon波束合成算法為核心，利用超聲面陣換能器接收腦部聚焦超聲治療過(guò)程中的空化信號(hào)，一方面可以彌補(bǔ)磁共振監(jiān)控和傳統(tǒng)主動(dòng)超聲成像方法的缺陷，另一方面可以解決現(xiàn)有無(wú)源成像方法在成像分辨率受限、成像偽影嚴(yán)重和三維空間信息丟失等方面的難題；通過(guò)改變聚焦超聲治療參數(shù)可以得到任意治療模式下的空化圖像，改變信號(hào)采集觸發(fā)參數(shù)可以得到任意幀頻的時(shí)間序列空化圖像，二者結(jié)合起來(lái)可反應(yīng)治療程度和治療效果，為臨床提供了一種腦部聚焦超聲治療過(guò)程中空化效應(yīng)的實(shí)時(shí)監(jiān)控手段，使得腦部治療的實(shí)時(shí)反饋與控制成為可能。

本發(fā)明具有以下優(yōu)點(diǎn)：

(1)針對(duì)磁共振成像與傳統(tǒng)主動(dòng)超聲成像無(wú)法與腦部聚焦超聲治療過(guò)程同步的問(wèn)題，通過(guò)設(shè)置面陣陣元變跡函數(shù)使其工作在不發(fā)射只接收模式，從而實(shí)現(xiàn)了腦部聚焦超聲治療中空化效應(yīng)的實(shí)時(shí)監(jiān)控，為腦部聚焦超聲治療的臨床應(yīng)用提供了一種實(shí)時(shí)反饋的手段。

(2)利用顱骨ct掃描數(shù)據(jù)對(duì)面陣換能器采集到的空化信號(hào)進(jìn)行校準(zhǔn)，在此基礎(chǔ)上對(duì)空化信號(hào)做沿著面陣兩個(gè)方向的魯棒的capon波束合成，可以克服顱骨遮擋引起的信號(hào)失真問(wèn)題，同時(shí)可以減弱來(lái)自于其他方向的干擾信號(hào)和空化微泡間的相互作用，使得三維無(wú)源成像的偽影被大幅度抑制。

(3)根據(jù)面陣陣元采集到的空化信號(hào)之間的相位差異引入了相位相干系數(shù)對(duì)波束合成算法進(jìn)行修正，在抑制了成像偽影的同時(shí)可以提高成像分辨率，使得對(duì)空化的三維空間定位更加精準(zhǔn)。

(4)高分辨及低偽影的三維無(wú)源成像方法使得對(duì)臨床治療的實(shí)時(shí)反饋更加直觀(guān)，更有利于對(duì)腦部治療過(guò)程的控制；提出的三維無(wú)源成像方法不僅可以用于腦部聚焦超聲治療的實(shí)時(shí)監(jiān)控，也可用于其他部位(如肝臟、腎臟等)和其他治療方式(如激光、微波等)的實(shí)時(shí)監(jiān)控。