超聲平面波經(jīng)顱成像相位校正方法?

宋亞龍 蘇暢,3 李倩巖 林偉軍,3

(1 中國科學(xué)院聲學(xué)研究所 北京 100190)

(2 中國科學(xué)院大學(xué) 北京 100049)

(3 北京市海洋深部鉆探測量工程技術(shù)研究中心 北京 100190)

0 引言

利用超聲穿透顱骨對腦組織進(jìn)行成像，在腦疾病臨床診斷和科研上具有重要應(yīng)用前景。然而當(dāng)超聲穿過顱骨時，由于顱骨的聲學(xué)特性(如密度、聲速、衰減)與周圍其他組織差異巨大[1]，導(dǎo)致超聲發(fā)生嚴(yán)重的衰減和畸變，給超聲經(jīng)顱腦成像帶來很大困難。由于超聲很難穿透顱骨，目前超聲經(jīng)顱成像一般透過顳骨、枕骨、下頜、眼眶、新生兒囟門等聲窗進(jìn)行，如經(jīng)顱彩超(Transcranial color-coded duplex sonography,TCCD)和經(jīng)顱多普勒(Transcranial Doppler,TCD)[2]。而近年來新型匹配材料和聲學(xué)超材料的發(fā)展[3]，可以增強(qiáng)超聲穿透顱骨的能量，使超聲穿透顱骨其他部位進(jìn)行治療或成像成為可能。此時由于顱骨帶來的超聲相位畸變會使超聲圖像出現(xiàn)分辨率下降、偽像等問題，需要加以校正以獲得準(zhǔn)確的顱內(nèi)圖像。

為了補(bǔ)償和校正顱骨對超聲相位的影響，Clement等[4?5]提出利用計(jì)算機(jī)斷層掃描/磁共振成像(CT/MRI)技術(shù)構(gòu)建非均勻顱骨模型，并通過理論建模和數(shù)值模擬計(jì)算經(jīng)顱聚焦超聲中陣列所需的相位補(bǔ)償。在應(yīng)用于腦疾病治療和神經(jīng)調(diào)控的經(jīng)顱聚焦超聲研究中，F(xiàn)ink[6]、Wu等[7]、Thomas等[8]提出的時間反轉(zhuǎn)法是公認(rèn)的最常用的解決聚焦偏差的方法。Aubry等[9]提出了結(jié)合CT圖像的虛擬點(diǎn)源時間反轉(zhuǎn)法，基于CT圖像獲取顱骨聲速模型，利用時域有限差分法對黏性流體介質(zhì)中的線性聲波方程[10]進(jìn)行數(shù)值求解獲得超聲相控陣各個陣元所接收到的信號，并使用時間反轉(zhuǎn)法計(jì)算用于補(bǔ)償?shù)南辔唬瑥亩鴮?shí)現(xiàn)超聲經(jīng)過顱骨在顱內(nèi)精準(zhǔn)聚焦。上述研究主要集中于高強(qiáng)度超聲聚焦(High intensity focused ultrasound,HIFU)治療和超聲神經(jīng)調(diào)控方面，對于腦部經(jīng)顱成像的研究相對較少。Ryan等[11]將類似的相位補(bǔ)償方法應(yīng)用于經(jīng)顱被動聲成像的接收波束形成，并對比了射線聲學(xué)理論和時間反轉(zhuǎn)數(shù)值模擬兩種補(bǔ)償方法對成像效果的影響。Sukhoruchkin等[12]從超聲回波信號中提取顱骨形態(tài)信息，用于經(jīng)顱超聲聚焦掃描成像中的相位補(bǔ)償。Lindsey等[13]基于Snell定律和分層介質(zhì)模型對三維經(jīng)顱成像進(jìn)行了時延校正。Chen等[14]利用射線聲學(xué)理論與理想合成孔徑聚焦技術(shù)(Ideal-Synthetic aperture focusing technique,I-SAFT)相結(jié)合，補(bǔ)償經(jīng)顱成像中的超聲相位，校正了圖像位置偏差并改善了圖像對比度。

另一方面，相比于傳統(tǒng)的聚焦掃描成像，平面波成像是一種高速成像方法，具有非常高的成像幀頻，有望滿足腦功能成像、實(shí)時三維成像、實(shí)時多普勒血流成像對高幀頻的需求。Macé等[15]利用平面波探測鼠的腦部，獲得了腦血流動態(tài)分布圖像。2019年，Du等[16]利用超聲發(fā)散波發(fā)射和深度學(xué)習(xí)技術(shù)進(jìn)行了經(jīng)顱成像，通過較小聲窗快速對顱腦內(nèi)的較大區(qū)域進(jìn)行組織結(jié)構(gòu)成像，并進(jìn)行了相關(guān)實(shí)驗(yàn)。這些工作中采用的顱骨較薄，因此并沒有對顱骨相位畸變進(jìn)行補(bǔ)償。胡陳文寶[17]在平面波經(jīng)顱高精度腦血流成像方面開展了動物實(shí)驗(yàn)，并對超聲測量顱骨厚度進(jìn)行了探索。

超聲平面波經(jīng)顱成像中超聲的發(fā)射方式和傳播路徑與超聲經(jīng)顱聚焦或被動聲成像有較大區(qū)別，因此其相位補(bǔ)償方法也有所不同。本文針對超聲平面波經(jīng)顱成像，利用已知的由CT圖像獲取的顱骨聲速模型，分別采用基于射線聲學(xué)近似的理論方法以及基于時間反轉(zhuǎn)的全波數(shù)值模擬方法，校正平面波經(jīng)顱入射和反向散射波穿過顱骨時的相位畸變，并利用數(shù)值仿真數(shù)據(jù)對比成像效果。

1 超聲平面波成像中的相干復(fù)合成像和經(jīng)顱相位畸變

超聲平面波成像時，利用換能器線陣同相位發(fā)射超聲脈沖，所形成的波陣面近似平面波，波前到達(dá)同一深度不同橫向位置的時間相同。平面波在聲阻抗發(fā)生變化的位置發(fā)生散射，此時散射點(diǎn)可看作被動點(diǎn)聲源，散射波以球面波形式傳播并被換能器接收，同一散射點(diǎn)回波到達(dá)換能器各陣元的時間與其之間距離成正比。一般的平面波成像假設(shè)介質(zhì)聲速均勻，根據(jù)成像區(qū)域各像素點(diǎn)與換能器陣元的距離關(guān)系計(jì)算聲傳播時間，對換能器陣列接收信號進(jìn)行延遲疊加(接收波束形成)，從而形成圖像。與傳統(tǒng)的聚焦掃描成像方法相比，平面波一次發(fā)射即可獲得完整的超聲圖像，大大提高了成像幀頻。單次發(fā)射的平面波成像分辨率和對比度較低，調(diào)整換能器各陣元發(fā)射相位，使換能器以不同角度發(fā)射平面波，并使用相干復(fù)合成像方法[18]，可以提高平面波成像的分辨率和對比度。

單次平面波成像的發(fā)射接收示意圖如圖1所示，換能器線陣沿x軸排列，各陣元以一定延遲發(fā)射，產(chǎn)生的平面波與x軸存在夾角α。以α小于π/2為例，并且以處于坐標(biāo)原點(diǎn)的陣元x1發(fā)射脈沖時為時間起始點(diǎn)，則陣元xk接收到散射點(diǎn)(x,y)產(chǎn)生的散射波的時間為

假設(shè)發(fā)射的平面波次數(shù)為M，與x軸的夾角分別為α1、α2、α3···αM，設(shè)RF(xk,t,αi)為不同平面波發(fā)射角度下接收陣元接收到的散射回波信號，將公式(1)代到RF(xk,t,αi)中，并且所有接收陣元進(jìn)行疊加，可以得到成像區(qū)域內(nèi)的像素點(diǎn)(x,y)處的像素值為

當(dāng)α大于π/2時，結(jié)果與α小于π/2時情況類似，不再過多闡述。

當(dāng)顱骨存在時，由于顱骨本身厚度和聲速分布不均勻，超聲沿不同路徑穿過顱骨不同位置時傳播時間不同，顱骨與周圍介質(zhì)聲速不同也使從不同角度入射并穿過顱骨的聲波傳播時間出現(xiàn)差異，因此超聲相位會發(fā)生畸變。入射平面波經(jīng)過顱骨后波陣面不再為平面，到達(dá)同一深度不同橫向位置處的波形相位有差別。而散射波經(jīng)過顱骨后也不再是球面波，同一散射目標(biāo)回波到達(dá)各陣元的時間與其間的直線距離不再為正比關(guān)系。

如果不做相位校正，成像結(jié)果會偏離實(shí)際散射點(diǎn)位置，并且出現(xiàn)偽像和分辨率下降等問題。因此，本文利用CT圖像獲取顱骨結(jié)構(gòu)和聲學(xué)參數(shù)模型，并采用基于近似射線聲學(xué)的理論方法和基于時間反轉(zhuǎn)的全波數(shù)值模擬方法，分別計(jì)算用于校正的相位偏差。

2 基于近似射線聲學(xué)的相位校正

不同于圖1所示的平面波發(fā)射接收，圖2顯示了平面聲波波前到達(dá)散射點(diǎn)和散射回波被換能器所接收這兩個過程中，聲波分別經(jīng)過了一段顱骨，也正是由于該兩段顱骨的存在導(dǎo)致成像結(jié)果出現(xiàn)偏差。事實(shí)上超聲在顱骨內(nèi)外表面會發(fā)生折射，由于顱骨聲速與周圍軟組織差異較大，超聲傳播路徑不應(yīng)是直線。在本文模型中，顱骨聲速分布為非均勻的，超聲在其中傳播路徑較復(fù)雜，因此在射線理論相位校正中，本文忽略了超聲折射造成的聲波傳播方向變化，簡單地用直線代替聲傳播路徑，并計(jì)算該路徑上由于顱骨聲速變化帶來的傳播時間差異。

圖1 平面波偏轉(zhuǎn)角度發(fā)射接收示意圖Fig.1 Schematic diagram of transmitting and receiving with plane wave of def lection angle

圖2 基于近似射線聲學(xué)理論的聲傳播路徑示意圖Fig.2 Schematic diagram of approximated acoustic propagation path based on ray theory

以圖2中成像區(qū)域內(nèi)像素點(diǎn)(x,y)為例，平面聲波由換能器發(fā)射，經(jīng)過顱骨傳播到像素點(diǎn)(x,y)時產(chǎn)生散射回波，再次經(jīng)過顱骨被換能器陣列中陣元mi所接收到，聲波從發(fā)射到接收過程中所經(jīng)過的兩次顱骨長度分別記為L1和L2。則兩段顱骨造成的時延偏差分別為

其中，cskull(l)代表的是顱骨的聲速，隨著顱骨中不同的位置而變化，cwater代表的是使用水的聲速代替成像時所使用的顱內(nèi)腦實(shí)質(zhì)的平均聲速(兩者數(shù)值接近)?？梢缘玫郊m正后的像素點(diǎn)(x,y)在陣元mi上的時延為

其中，H1和H2分別代表聲波從換能器陣列發(fā)射到像素點(diǎn)(x,y)的傳播路徑長度和散射回波從像素點(diǎn)(x,y)傳播到陣元mi的傳播路徑長度。將公式(5)帶入到公式(2)中就可得到經(jīng)過相位補(bǔ)償后的像素點(diǎn)(x,y)處的像素值。使用近似射線法進(jìn)行相位補(bǔ)償比較直觀，易于操作，但本文沒有考慮實(shí)際聲傳播路徑上的折射效應(yīng)，會導(dǎo)致一定誤差。

3 基于時間反轉(zhuǎn)和數(shù)值計(jì)算的相位校正和成像

虛擬點(diǎn)源時間反轉(zhuǎn)法在進(jìn)行超聲經(jīng)顱聚焦過程中有著重要的作用，利用時間反轉(zhuǎn)計(jì)算顱骨導(dǎo)致的相位偏差，并在超聲發(fā)射時加以補(bǔ)償，可以使超聲精準(zhǔn)地聚焦至顱內(nèi)目標(biāo)處。同樣地，在平面波經(jīng)顱腦成像中，也可以利用時間反轉(zhuǎn)方法來計(jì)算和補(bǔ)償顱骨造成的聲波信號相位畸變，在超聲平面波入射至顱內(nèi)和散射波經(jīng)顱骨向外傳播兩個過程中，采用時間反轉(zhuǎn)方法分別計(jì)算所需的相位補(bǔ)償量，并分別在發(fā)射和接收過程中進(jìn)行補(bǔ)償處理。

3.1 時間反轉(zhuǎn)平面波發(fā)射相位校正

利用基于虛擬線陣的時間反轉(zhuǎn)數(shù)值模擬，可以計(jì)算平面波經(jīng)顱傳播時由顱骨導(dǎo)致的相位畸變，進(jìn)而在超聲發(fā)射時調(diào)控各陣元的發(fā)射相位，使超聲波陣面在進(jìn)入顱骨后以類似平面波的形式傳播。發(fā)射相位調(diào)控的過程與時間反轉(zhuǎn)超聲經(jīng)顱聚焦類似，區(qū)別在于所用虛擬聲源不是點(diǎn)聲源，而是一組平面線陣。如圖3(a)所示，從顱骨內(nèi)的虛擬線陣發(fā)射聲波，將顱外實(shí)際用于成像的換能器陣作為接收，利用數(shù)值仿真計(jì)算顱外換能器各陣元所在位置的聲波時域信號，記作fn(t)，其中下標(biāo)n為陣元編號。從這些信號提取相位差，用于發(fā)射校正，具體操作為：設(shè)置某一陣元接收到的聲波時域信號為參考信號，記為rref(t)，將其他陣元接收到的時域信號與參考信號進(jìn)行互相關(guān)處理[19]，并對參考信號進(jìn)行自相關(guān)處理，從而可得到

其中，Rref,n表示參考信號rref(t)和第n個陣元接收到的時域聲場信號rn(t)之間的互相關(guān)函數(shù)，Rref,ref是參考信號的自相關(guān)函數(shù)，τref,n是互相關(guān)信號的時間延遲，τref是自相關(guān)函數(shù)的時間延遲。找出分別使Rref,n(τref,n)和Rref,ref(τref)取最大值的τref,n和τref，就可以求得每個陣元的初始時延補(bǔ)償。

換能器各陣元利用以上時延補(bǔ)償調(diào)控相位并發(fā)射，如圖3(b)所示，則聲波穿過顱骨時的相位偏差得以補(bǔ)償，在顱內(nèi)以近似平面波的形式傳播，并且可以準(zhǔn)確地獲得用以成像的時間延遲。圖4給出了各陣元同相位發(fā)射時(無發(fā)射角度偏轉(zhuǎn))分別使用近似射線法和全波數(shù)值模擬計(jì)算得到的由虛擬陣列至實(shí)際陣列的各陣元時延偏差，以及它們之間的差值，可以看出，在發(fā)射端進(jìn)行相位補(bǔ)償時，兩種方法之間的相位補(bǔ)償差值不超過0.2μs。由于換能器孔徑相對較小，顱骨表面弧度較小，無角度偏轉(zhuǎn)的平面波可看作近似于垂直入射，因此用直線代替聲傳播路徑的近似射線法在平面波入射段誤差較小。

圖3 虛擬聲源時間反轉(zhuǎn)法用于平面波經(jīng)顱成像示意圖Fig.3 Schematic diagram of plane wave transcranial imaging using time reversal method based on virtual source

圖4 使用近似射線法和時間反轉(zhuǎn)法的時延偏差及其差值Fig.4 Time delay deviation and its difference between ray method and time reversal method

將時間反轉(zhuǎn)法得到的時延偏差用于發(fā)射相位校正，利用數(shù)值模擬仿真驗(yàn)證效果，圖5給出了不加相位校正和時間反轉(zhuǎn)相位校正發(fā)射后，顱內(nèi)距換能器18 mm處接收到的時域聲波信號，可見當(dāng)聲波穿過顱骨后，校正后的超聲波前可近似看作平面波。

圖5 距離換能器18 mm處的聲波接收信號Fig.5 Received acoustic signals at 18 mm from transducer

3.2 時間反轉(zhuǎn)接收相位校正

平面波經(jīng)顱成像接收端的相位校正相對于發(fā)射端來講復(fù)雜得多，需要對成像區(qū)域內(nèi)的每個像素點(diǎn)進(jìn)行相位補(bǔ)償，這就造成了巨量的計(jì)算時間和計(jì)算資源，特別是成像區(qū)域較大時，不加處理使用傳統(tǒng)的時間反轉(zhuǎn)法更是不符合實(shí)際需求的。Ryan等[11]將被動聲成像焦點(diǎn)處的相位補(bǔ)償值用于整個成像區(qū)域，這是在成像區(qū)域集中于焦點(diǎn)附近時采用的近似處理。但平面波聲傳播路徑不同，且成像區(qū)域相對較大，必須對不同位置分別計(jì)算相位補(bǔ)償量。為了減小時間反轉(zhuǎn)的數(shù)值計(jì)算量，本文提出了基于虛擬線陣的理論-數(shù)值混合算法。由于顱內(nèi)的腦實(shí)質(zhì)部分可以近似為均勻介質(zhì)，超聲在顱腦以內(nèi)的傳播過程可以用解析方法計(jì)算，而超聲穿透顱骨傳播的過程則由數(shù)值方法求解[20]。為此，在顱腦內(nèi)靠近顱骨的位置放置一組虛擬線陣，如圖6所示。首先，把顱內(nèi)一個像素點(diǎn)當(dāng)作被動點(diǎn)聲源，點(diǎn)聲源發(fā)射的聲波以球面波形式傳播至虛擬線陣，每個虛擬陣元處的聲波幅度和相位不同，此過程可解析求解。接下來，虛擬線陣各陣元以此相位延遲和幅度發(fā)射發(fā)散波，穿過顱骨到達(dá)顱外的換能器陣列，可以近似地代替點(diǎn)聲源產(chǎn)生的聲場，此過程可通過數(shù)值模擬結(jié)合聲場疊加原理來求解。

圖6(a)中假設(shè)顱腦內(nèi)側(cè)的虛擬陣列B陣元為Nb個，成像區(qū)域內(nèi)任一像素點(diǎn)S可看作一個被動點(diǎn)聲源，該點(diǎn)聲源S輻射的球面波到達(dá)虛擬陣列B，陣列B各陣元接收到的聲波信號為Ais(t?τi)，其中i為陣元序號，Ai和τi可由解析解獲得。圖6(b)中虛擬陣列B各陣元以Ais(t?τi)為發(fā)射信號進(jìn)行發(fā)射，顱骨外的陣列A接收的信號，應(yīng)與點(diǎn)聲源S直接發(fā)射時近似相等。實(shí)際操作中，陣列B中每個陣元i單獨(dú)發(fā)射，陣列A陣元接收到信號Qj(i,t)，j為陣元序號，根據(jù)聲場疊加原理，陣列A中陣元j接收的聲源S發(fā)射的聲波信號rcvj(t)可由Qj(i,t)通過延遲加權(quán)疊加來獲?。?/p>

當(dāng)介質(zhì)中無顱骨存在時，陣列A各陣元接收的聲源S發(fā)射的聲波信號refj(t)作為參考信號，與rcvj(t)做互相關(guān)處理，得到的時間延遲即接收波束形成中的顱骨延時補(bǔ)償量，可用于成像中的相位補(bǔ)償。

按傳統(tǒng)時間反轉(zhuǎn)方法，對成像區(qū)域內(nèi)每個像素點(diǎn)(x,y)都需要數(shù)值計(jì)算聲傳播過程，以獲得各陣元的延遲補(bǔ)償量。而利用以上的虛擬線陣方法，只需進(jìn)行Nb次聲傳播的數(shù)值計(jì)算，大大減少了計(jì)算量，同時由于虛擬陣列貼近顱骨，聲傳播的數(shù)值計(jì)算區(qū)域也可顯著減小。

圖6 接收相位校正示意圖Fig.6 Receiving phase correction diagram

4 數(shù)值仿真

根據(jù)Aubry等[9]建立的顱骨的聲參數(shù)與顱骨亨氏值之間的關(guān)系可以得到顱骨的聲參數(shù)：

其中，HU是顱骨亨氏值，可以由顱骨的CT文件中獲得；φ是指顱骨的孔隙率；ρmin和ρmax是顱骨中的最小密度和最大密度；cmin和cmax是顱骨中的最小聲速和最大聲速；本文選取了顱骨中較厚的一段用來構(gòu)建計(jì)算模型，如圖7所示。仿真中使用的顱骨和腦實(shí)質(zhì)的具體聲參數(shù)如表1所示。

表1 仿真使用的聲參數(shù)Table 1 Acoustic parameters used in the simulation

圖7 基于顱骨CT文件的顱骨模型Fig.7 Skull model based on skull CT file

本文根據(jù)獲得的顱骨非均勻參數(shù)模型建立整體的計(jì)算模型如圖8所示，仿真中所使用的是平面線陣，陣列放置在顱骨外距顱骨大約2 mm處，發(fā)射中心頻率為2.4 MHz，陣元個數(shù)為64個，孔徑為25.2 mm，陣元間距約為0.4 mm(略微大于半波長)，陣列發(fā)射的聲波形式為

陣列與顱骨之間使用水進(jìn)行耦合。在顱腦內(nèi)近似均勻的腦實(shí)質(zhì)環(huán)境內(nèi)放置8個散射目標(biāo)點(diǎn)，目標(biāo)點(diǎn)橫向間距2 mm，縱向間距5 mm，分布在距換能器20～35 mm范圍內(nèi)，目標(biāo)點(diǎn)的聲速和密度略大于腦實(shí)質(zhì)的聲速和密度。

圖8 平面波經(jīng)顱成像計(jì)算模型Fig.8 Calculation model of plane wave transcranial imaging

在進(jìn)行仿真計(jì)算時，本文基于流體內(nèi)的線性聲波方程，使用fortran語言編寫了空間四階精度、時間二階精度的二維交錯網(wǎng)格時域有限差分法程序，計(jì)算中沒有考慮顱骨中的橫波，并且忽略了介質(zhì)衰減。計(jì)算區(qū)域?yàn)閤軸方向(橫向方向)40 mm、y軸方向(深度方向)60 mm。計(jì)算空間步長為0.08 mm，約為水中波長的1/10，時間步長為0.01μs，滿足Courant-Friedrich計(jì)算穩(wěn)定性條件。

利用仿真數(shù)據(jù)，對x軸方向(橫向方向)?10～10 mm、y軸方向(深度方向)15～40 mm的區(qū)域，用前文所述的平面波成像方法和兩種相位補(bǔ)償方法進(jìn)行了成像，圖9給出了單角度平面波發(fā)射時(聲傳播方向垂直于發(fā)射陣列)的成像結(jié)果。其中，圖9(a)為無顱骨存在時的成像結(jié)果，圖9(b)為有顱骨但未做相位補(bǔ)償?shù)膱D像，圖9(c)、圖9(d)分別為近似射線法和時間反轉(zhuǎn)法補(bǔ)償后的圖像。圖中藍(lán)色點(diǎn)代表散射目標(biāo)點(diǎn)原本所處位置，圖9(d)中的黃色陣列代表時間反轉(zhuǎn)補(bǔ)償時顱骨內(nèi)虛擬陣列放置位置(距發(fā)射陣列18 mm)。此外，單角度平面波成像往往帶來比較明顯的偽像，有顱骨存在時進(jìn)一步增加了偽像，如圖9(b)、圖9(c)、圖9(d)中紅色虛線圈出區(qū)域所示(部分偽像)。

圖9所顯示的結(jié)果是單角度平面波發(fā)射時的成像結(jié)果，不可避免地會有平面波成像分辨率差、對比度差等缺陷，實(shí)際平面波成像的應(yīng)用中一般采用多角度平面波發(fā)射和相干復(fù)合法提高圖像質(zhì)量[18]。如本文第1小節(jié)所述，從?10?～+10?、間隔1?發(fā)射不同角度平面波，共發(fā)射21次。對每個角度發(fā)射后接收的超聲信號分別使用兩種方法進(jìn)行相位補(bǔ)償和成像，最后將不同發(fā)射角度的結(jié)果進(jìn)行相干復(fù)合疊加，得到最終的成像圖，如圖10和圖11所示。圖10和圖11分別給出了用近似射線法和時間反轉(zhuǎn)法進(jìn)行相位補(bǔ)償后，不同相干復(fù)合次數(shù)下成像的結(jié)果。對不同復(fù)合次數(shù)，平面波發(fā)射角度間隔不變，最大偏轉(zhuǎn)角度不同，圖10(a)、圖11(a)為5次平面波發(fā)射相干復(fù)合的成像結(jié)果，圖10(b)、圖11(b)為11次相干復(fù)合的圖像，圖10(c)、圖11(c)為21次平面波發(fā)射相干復(fù)合后的圖像。

圖9 單角度平面波發(fā)射時成像結(jié)果Fig.9 Imaging results of single angle plane wave emission

圖10 近似射線法補(bǔ)償后不同復(fù)合次數(shù)的成像結(jié)果Fig.10 Imaging results of different complex times after phase compensation with approximated ray method

圖11 時間反轉(zhuǎn)法補(bǔ)償后不同復(fù)合次數(shù)的成像結(jié)果Fig.11 Imaging results of different complex times after phase compensation with time reversal method

5 討論

對比單角度平面波發(fā)射時的圖像，從圖9(a)、圖9(b)可以直觀地看出，相對于無顱骨存在時，在不進(jìn)行相位補(bǔ)償?shù)那闆r下，顱骨的存在使得圖像產(chǎn)生了嚴(yán)重的偏差。首先，像點(diǎn)在深度方向上偏離實(shí)際位置，向靠近換能器方向移動，這主要是顱骨聲速遠(yuǎn)大于周圍介質(zhì)聲速造成的。其次，圖像分辨率下降，在水平方向上相對更為顯著，間距為2 mm的兩個散射點(diǎn)完全不能分辨，在圖中幾乎顯示為連續(xù)界面。此外，圖像對比度也有明顯下降。

根據(jù)圖9(c)、圖9(d)，使用兩種相位補(bǔ)償方法都能夠顯著改善成像質(zhì)量，表2給出了量化的成像質(zhì)量對比，包括8個目標(biāo)散射點(diǎn)圖像的位置偏差(?L)、像點(diǎn)水平方向上幅度下降至?3 dB的寬度(W)以及像點(diǎn)最亮處相對于其兩側(cè)偽像幅度的對比度(C)。首先，兩種相位校正方法都起到了校正像點(diǎn)位置偏差的作用：當(dāng)沒有進(jìn)行相位補(bǔ)償時，圖像中亮斑與實(shí)際散射點(diǎn)的位置偏差平均約為3.8 mm；使用近似射線法進(jìn)行相位校正之后，像點(diǎn)位置平均偏差為0.4 mm，減小了90%；而在使用時間反轉(zhuǎn)法進(jìn)行相位校正之后，位置偏差幾乎可以忽略不計(jì)，所成出的圖像可以準(zhǔn)確地反映散射目標(biāo)點(diǎn)所處位置。其次，相位校正對圖像橫向分辨率有明顯改善作用：未進(jìn)行相位補(bǔ)償操作時，像點(diǎn)亮斑橫向?3 dB的平均寬度約為1.7 mm，遠(yuǎn)遠(yuǎn)大于模型中設(shè)置的散射目標(biāo)點(diǎn)的直徑0.64 mm，約為水中超聲波長的3倍；分別使用近似射線法和時間反轉(zhuǎn)法進(jìn)行相位補(bǔ)償處理后?3 dB的平均寬度為1.3 mm和0.9 mm。最后，補(bǔ)償后的圖像對比度也得到顯著提高，近似射線法和時間反轉(zhuǎn)法得到的像斑與偽像幅度之比分別提高至未校正時的1.75倍和2.05倍。綜上所述，兩種相位方法都可以有效地減小顱骨造成的像的位置偏差，提高成像分辨率和對比度。使用時間反轉(zhuǎn)法進(jìn)行相位校正的精度和效果好于近似射線法，與經(jīng)顱聚焦超聲[9]、經(jīng)顱被動聲成像[11]等研究中結(jié)果一致。此外，本文使用的近似射線法時并沒有考慮非均勻顱骨造成的聲波折射，與嚴(yán)格的射線聲學(xué)方法[14]相比可能帶來了更多誤差。

對比圖10和圖9(c)、圖11和圖9(d)，不難看出，在經(jīng)顱成像過程中，平面波相干復(fù)合成像可以有效抑制偽像的產(chǎn)生，提高圖像分辨率和對比度。此外，同等復(fù)合次數(shù)下，使用時間反轉(zhuǎn)法相位補(bǔ)償后的平面波經(jīng)顱復(fù)合成像情況明顯好于使用近似射線法相位補(bǔ)償后的平面波經(jīng)顱復(fù)合成像情況：前者需要復(fù)合20次才能獲得比較不錯的成像結(jié)果，而后者復(fù)合10次，甚至是5次就可以獲得理想的成像效果。

表2 成像結(jié)果對比Table 2 Comparison of imaging results

雖然使用時間反轉(zhuǎn)法的相位補(bǔ)償效果要好于近似射線法的補(bǔ)償效果，但是其計(jì)算過程更為復(fù)雜，所需的計(jì)算時間遠(yuǎn)遠(yuǎn)超過近似射線法。時間反轉(zhuǎn)法需要在成像發(fā)射時調(diào)整各陣元發(fā)射相位，對成像設(shè)備及其控制提出了更高要求，在接收相位校正時對虛擬陣列每個陣元做全波聲傳播過程的數(shù)值模擬計(jì)算，又需要大量計(jì)算資源和時間。以本研究中所用仿真數(shù)據(jù)為例，對x軸方向(橫向方向)?10～10 mm、y軸方向(深度方向)15～40 mm的區(qū)域成像，像素點(diǎn)數(shù)為400×500。用近似射線法進(jìn)行相位補(bǔ)償時，讀取顱骨模型和仿真數(shù)據(jù)、相位校正和成像過程總用時約45 s。而利用時間反轉(zhuǎn)法對相同區(qū)域成像時，采用64陣元的虛擬線陣，即需要進(jìn)行64次二維聲場模擬，盡管模擬計(jì)算中將虛擬線陣靠近顱骨從而盡量縮小了聲場計(jì)算區(qū)域，并使用了4個進(jìn)程并行計(jì)算以加快速度，一次成像的總用時仍需約3 h，且在計(jì)算過程中需要約1 GB的內(nèi)存來存儲所需的雙精度參數(shù)。以上模擬和成像計(jì)算均采用自行編寫的fortran軟件，在配有4核主頻3.10 GHz處理器的PC機(jī)上進(jìn)行。由于時間反轉(zhuǎn)法計(jì)算時間過長，暫時還難以在實(shí)際的顱腦成像中應(yīng)用，并且需要高性能計(jì)算機(jī)和快速模擬算法的支持。此外，兩種相位校正方法都依賴于準(zhǔn)確的顱骨聲速先驗(yàn)?zāi)Ｐ停@在實(shí)際成像中是難以獲取的。因此研究加快顱腦聲場模擬計(jì)算速度的算法、尋求不依賴顱骨先驗(yàn)?zāi)Ｐ蛯?shí)時獲取顱骨聲速的方法或用簡化顱骨模型提高顱腦成像精度，是使超聲經(jīng)顱成像走向?qū)嵱玫膸讉€可能的發(fā)展方向。

6 結(jié)論

顱骨的存在使得在平面波超聲經(jīng)顱成像中的超聲傳播產(chǎn)生了嚴(yán)重的相位畸變，從而導(dǎo)致成像結(jié)果存在位置偏差、分辨率和對比度低等問題。本文分別使用近似射線法和基于虛擬聲源的時間反轉(zhuǎn)法對顱骨造成的相位畸變進(jìn)行了補(bǔ)償，并利用交錯網(wǎng)格時域有限差分求解無黏線性聲波方程對穿顱聲場進(jìn)行了建模仿真和成像驗(yàn)證。結(jié)果表明：(1)相對于無補(bǔ)償時的成像結(jié)果，無論使用近似射線法還是時間反轉(zhuǎn)法，都能夠有效地校正因顱骨造成的相位畸變，從而減小像點(diǎn)位置偏差、分辨率、對比度低等問題。(2)使用時間反轉(zhuǎn)和近似射線法的仿真結(jié)果存在微小的偏差，主要是由于在使用射線法時并沒有考慮顱骨造成的聲波折射，從而在確定聲傳播路徑時存在一定誤差。(3)時間反轉(zhuǎn)法成像的精度要好于近似射線法，但所需的計(jì)算資源和時間都要遠(yuǎn)遠(yuǎn)大于近似射線法。(4)無論近似射線法還是時間反轉(zhuǎn)法，經(jīng)過平面波多角度發(fā)射和相干復(fù)合處理后都能夠一定程度上提高成像的對比度和分辨率。

致謝感謝浙江大學(xué)醫(yī)學(xué)院附屬第四醫(yī)院提供的顱骨CT掃描文件。