超聲平面波經顱成像相位校正方法?

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

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

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

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

0 引言

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

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

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

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

1 超聲平面波成像中的相干復合成像和經顱相位畸變

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

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

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

當α大于π/2時，結果與α小于π/2時情況類似，不再過多闡述。

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

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

2 基于近似射線聲學的相位校正

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

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

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

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

其中，cskull(l)代表的是顱骨的聲速，隨著顱骨中不同的位置而變化，cwater代表的是使用水的聲速代替成像時所使用的顱內腦實質的平均聲速(兩者數值接近)。可以得到糾正后的像素點(x,y)在陣元mi上的時延為

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

3 基于時間反轉和數值計算的相位校正和成像

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

3.1 時間反轉平面波發射相位校正

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

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

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

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

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

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

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

3.2 時間反轉接收相位校正

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

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

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

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

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

4 數值仿真

根據Aubry等[9]建立的顱骨的聲參數與顱骨亨氏值之間的關系可以得到顱骨的聲參數：

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

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

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

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

陣列與顱骨之間使用水進行耦合。在顱腦內近似均勻的腦實質環境內放置8個散射目標點，目標點橫向間距2 mm，縱向間距5 mm，分布在距換能器20～35 mm范圍內，目標點的聲速和密度略大于腦實質的聲速和密度。

圖8 平面波經顱成像計算模型Fig.8 Calculation model of plane wave transcranial imaging

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

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

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

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

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

圖11 時間反轉法補償后不同復合次數的成像結果Fig.11 Imaging results of different complex times after phase compensation with time reversal method

5 討論

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

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

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

表2 成像結果對比Table 2 Comparison of imaging results

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

6 結論

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

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