動脈粥樣硬化及血栓的形成與動脈管腔內的血流動力學參數變化密切相關。然而,目前普遍應用的超聲多普勒成像技術不能精確測量復雜血流流場信息。本文提出了一種基于超聲造影微泡的超聲粒子成像測速技術,能夠獲得二維血管全流場的信息,且不依賴于聲束與速度向量之間的夾角。本研究使用10只健康的大鼠,應用超聲粒子成像測速法計算其左側頸動脈的全場血流速度分布,結合多項式擬合的方法計算腔內速度梯度變化,可評估血管在不同血流速度下受到的剪切力。實驗結果表明超聲粒子成像技術測速結果與多普勒測速結果吻合。利用四次多項式擬合的正常血管流速模型效果最佳。統計血管壁周期的平均剪切力結果與該位置平均血流速度成正相關。超聲粒子成像技術能夠無創、實時地評估全流場血流動力學變化情況,為進一步探討狹窄血管流場中復雜多變的剪切力分布異常奠定基礎。
引用本文： 朱懿恒, 錢明, 牛麗麗, 鄭海榮, 盧廣文. 基于超聲粒子成像技術的大鼠動脈模型血流剪切力測量方法. 生物醫學工程學雜志, 2014, 31(6): 1355-1360. doi: 10.7507/1001-5515.20140257 復制
引言
近年來心血管疾病的發病率和死亡率逐年提升，占人類非傳染性疾病總死因的第一位，嚴重威脅著人類健康[1-2]。動脈血管內血流動力學的改變會影響血管壁內皮細胞的形態和功能[2]。究其根本，是由于血液的流動會對血管壁產生包括切應力(剪切力)、周應力和壓應力等在內的各種作用力，從而對血管內皮細胞進行不斷的刺激。其中切應力被認為是與動脈粥樣硬化發生及斑塊破裂最為相關的力學因素[3]。測量血流剪切力要首先測量血管內血流流速分布，并通過剪切力與血流流速梯度之間的正比關系而求出。目前，測量血管血流速度常見的影像學方法是超聲多普勒和磁共振成像(magnetic resonance imaging,MRI)。超聲多普勒只能提供一維的速度信息。為了測量剪切力，需假設血管內流速是標準的層流分布且流速具備拋物線狀分布特征。而這種流速特征只有在形狀規則的血管內才能成立。對于復雜形狀(彎曲、分叉)的血管內的流速分布就會有較大的誤差。而且超聲多普勒技術對聲束與速度向量之間的夾角依賴性大，亦可能產生一定的誤差。MRI具有較高的空間分辨率,但時間分辨率較低，且檢查費用昂貴[4-5]。
超聲粒子成像測速是一種基于粒子成像測速法提出的實時、動態、無創、高分辨率的超聲測速方法。它可以得到精細的二維流場速度矢量、速度梯度分布，進而計算流體剪切力分布[5]。相關研究指出[6]，應用超聲粒子成像測速方法能夠精確計算血管仿體中流體流場的速度分布以及大鼠動脈和靜脈的血流速度。
本文應用超聲粒子成像測速技術獲取大鼠左側頸總動脈(left carotid artery,LCA)血流速度矢量分布，并結合多項式擬合方法計算血管剪切力分布。采集10只大鼠的周期性動脈血流信息，并將超聲粒子成像測速結果與多普勒測速結果進行對比。統計動脈血管不同位置的血流速度和剪切力周期性變化，對于研究血管壁內膜增生和斑塊生長、破裂的機制有重要的意義。
1 材料與方法
超聲粒子圖像測速技術是利用造影微泡對超聲波的強烈反射，并通過高幀頻超聲系統對造影區域采集圖像序列，然后對連續的兩幅圖像進行圖像相關處理，計算圖像中微泡的位移，進而獲得流場的速度和剪切力分布信息。如圖 1所示，超聲粒子圖像測速系統主要由動物超聲系統、圖像序列存儲器、圖像處理算法組成。
圖1
超聲粒子成像測速系統
Figure1.
Ultrasonic particle imaging velocimetry system
1.1 超聲系統
實驗采用由加拿大VisualSonics公司生產的Vevo2100超聲實時分子影像系統。與傳統的超聲系統相比，Vevo2100具備超高的空間分辨率(高達30 μm)和時間分辨率(大于1 000幀/s)。該系統采用了超高頻率的電子線陣探頭(探頭型號為MS400，頻率高達40 MHz)，適用于觀察大鼠微小血管和快速心跳頻率下的造影血流。
1.2 實驗動物及圖像數據采集
10只健康雄性SD大鼠,體重(340±10)g。大鼠被麻醉后，去除左側頸部毛發并固定在試驗臺，使用Vevo2100超聲實時分子影像系統(探頭頻率為40 MHz，幀頻為200 Hz)檢測大鼠左側頸總動脈。超聲檢查始終由同一醫師操作，檢查時超高頻探頭使用支架固定，測量多普勒血流速度。使用5 mL生理鹽水配制超聲造影劑(六氟化碳微泡，聲諾維，SonoVue)并充分搖勻。用1 mL注射器抽取0.2 mL超聲造影劑，經大鼠右側股靜脈注入，再用0.1 mL生理鹽水沖管。同時使用Vevo2100超聲實時分子影像系統對大鼠在造影劑注射時、注射5 min后、注射10 min后分別采集1 000幀圖像，如圖 2所示。
圖2
大鼠LCA的B超圖像
超聲造影劑
Figure2. B-mode images of the rat LCAUCAs: ultrasound contrast agents
1.3 超聲粒子成像測速算法
本文所應用的超聲粒子成像測速法把圖像互相關算法與多次迭代算法、亞像素邊緣檢測方法相結合，并通過濾波和插值的方法消除偽向量，從而改善傳統超聲粒子成像測速算法的測量精度和效率[7]。首先輸入連續的兩幅圖像，設置迭代次數K，再根據圖像互相關算法和亞像素峰值擬合方法計算位移矢量，經插值濾波后得到第一次迭代的結果。若迭代沒有完成，根據前一次迭代的結果進行詢問窗口變形，再重復上述迭代過程。迭代完成后減小詢問窗口大小，應用圖像互相關算法計算位移矢量，然后利用多項式擬合的方法求剪切力，算法流程如圖 3所示。
圖3
算法流程圖
Figure3.
Flow chart of the algorithm
1.3.1 圖像互相關算法
為了獲得流場速度分布，利用圖像互相關算法處理記錄的連續兩幅超聲粒子圖像。首先選擇合適的詢問窗口尺寸，然后超聲粒子圖像被分割成若干個小的詢問窗口，利用互相關函數來衡量兩個詢問窗口之間的灰度相似程度。根據相關文獻的方法[8]，互相關平面內峰值位置直接反映了詢問窗口的位移，再根據連續兩幅圖像的時間差和圖像放大倍數，可獲得該詢問窗口內的流體速度向量。對每一個詢問窗口重復以上步驟就可以獲得二維全流場速度矢量圖。以下是圖像互相關算法公式(1)：
| ${{R}_{ab}}(i,j)=\frac{\sum\limits_{x=0}^{M}{\sum\limits_{y=0}^{N}{[a(i,j)-][b(x-i,y-j)-]}}}{{{\sigma }_{a}}{{\sigma }_{b}}}~,$ |
其中M和N是詢問窗口的維數。a(i,j)和b(i,j)表示灰度值，和表示a和b的平均灰度值，σa和σb是標準差。但是圖像互相關方法在測量亞像素的位移以及粒子團在高速度梯度下的位移時誤差較大，因此本文使用亞像素邊緣檢測方法和多次迭代的方法彌補圖像互相關算法的不足。
1.3.2 亞像素邊緣檢測方法
圖像的灰度是以離散函數的像素單元組成，而位移向量也是以像素為單位。亞像素方法可以評估亞像素的位移，從而提高位移的測量精度。亞像素高斯峰值擬合方法是亞像素分析中常用的方法，以下是亞像素高斯峰值擬合的公式(2)：
| $\left\{ \begin{align} & \Delta x=i+\frac{ln(R(i-1,j))-ln(R(i+1,j))}{2ln(R(i-1,j))-4ln(R(i,j))+2ln(R(i+1,j))} \\ & \Delta y=j+\frac{ln(R(i,j-1))-ln(R(i,j+1))}{2ln(R(i,j-1))-4ln(R(i,j))+2ln(R(i,j+1))} \\ \end{align} \right.,$ |
其中R是互相關算法函數，(i,j)是互相關函數R取得最大值對應的坐標。
1.3.3 多次迭代算法
詢問窗口的粒子團在高速度梯度下流動會產生幾何形變，需要通過詢問窗口變形的方法解決。運行K次迭代算法，第K次迭代時詢問窗口的變形根據第K-1次迭代計算得的位移場Uk-1來獲得。而詢問窗口的變形是通過利用圖像像素點的位移來實現，對于互相關窗口的兩個圖像灰度函數a(i,j)與b(i,j)，利用算法重建后:
| $a\left( r \right)=ar-\frac{{{U}_{K-1}}}{2}-\frac{U{{\prime }_{k-1}}r}{2},$ |
| $b\left( r \right)=br+\frac{{{U}_{K-1}}}{2}-\frac{U{{\prime }_{k-1}}r}{2},$ |
其中r=(i,j),位移梯度U可通過中心有限差分獲得。對經過移位后的圖像再次利用互相關算法和亞像素高斯峰值擬合法計算函數R取得最大值所對應的坐標Rmax，而迭代后的位移UK為：
| ${{U}_{K}}={{U}_{K-1}}+{{R}_{max}}$ |
1.4 多項式擬合法計算血管剪切力
血液在血管腔內流動時，血液與血管壁接觸面上因相對運動產生力F，使血液產生速度變化du。由于液體的粘性將此力層層傳遞，各層液體也相應運動，形成速度梯度du/dr，稱剪切率。設血液與血管壁的接觸面積為A，F/A稱為剪切應力。剪切率與剪切應力具有如下關系：
| $\left( F/A \right)=\eta (du/dr)$ |
其中dr是層厚，η為液體黏度。血管壁受到的剪切力與血流的剪切率成正比，因此利用超聲粒子成像測速計算血流的速度梯度就可評估血管受到的剪切力大小。
血管剪切力是根據血管血流流速具備拋物線狀分布特征來計算的，先應用多項式擬合法擬合拋物線狀分布的血流流速分布，再根據擬合后結果計算剪切力。使用最小二乘法擬合多項式，假設給定數據點(xi,yi)(i=0,1,…,m)，Φ為所有次數不超過n(n≤m)的多項式構成的函數類，求，使得
| $I=\sum\limits_{i=0}^{m}{{{[{{p}_{n}}({{x}_{i}})-{{y}_{i}}]}^{2}}}=\sum\limits_{i=0}^{m}{\left( \sum\limits_{i=0}^{N}{{{a}_{k}}{{x}_{i}}^{k}-{{y}_{i}}^{2}} \right)}=min,$ |
再根據血流流速分布情況確定擬合多項式次數，滿足式(7)的pn(x)稱為最小二乘擬合多項式。
2 結果
本算法基于MATLAB R2009a調試運行，使用24×10至48×16像素大小的詢問窗，0.5~0.75的窗口重疊率，每只大鼠選取其血管微泡濃度適當的圖像序列(約200幀)進行處理。如圖 4所示,選取其中一只大鼠的造影圖像序列，圖 4(a)中黃色虛線區域是大鼠左側頸動脈,在圖像中截取感興趣區域(region of interest,ROI)，選擇像素大小為24×10的詢問窗和50%的窗口重疊率，計算該區域的血流速度分布。圖 4(b)中紅色箭頭的方向和長度分別表示大鼠頸動脈的血流速度方向和大小，并疊加在其對應的血管血流影像區域。如圖 4(c)所示，在血管中心區域血流速度較大，越靠近管壁兩邊血流速度越小。
圖4
應用超聲粒子成像測速技術的大鼠LCA二維血流速度圖
(a)管腔內注射UCAs的大鼠LCA的B超圖像，紅色方框是LCA所在的ROI；(b)ROI的速度矢量分布圖并疊加在其B超圖像上；(c)ROI血流速度梯度分布圖
Figure4. Two-dimensional flow velocity map of rat LCA calculated by ultrasonic particle image velocimetry(a) the B-mode images of the rat LCA with UCAs,and the red box marked the ROI of LCA; (b) velocity vectors map of ROI in the B-mode images; (c) flow velocity gradient map of ROI
為了驗證超聲粒子成像測速方法的可靠性，對連續三個心動周期(約200 ms)超聲粒子成像測速的最大速度取平均值，獲得一個周期的超聲粒子成像測速曲線，并與超聲多普勒測速對比。如圖 5所示，應用本文方法處理其中一只大鼠的數據并與多普勒測速結果對比，超聲粒子成像測速血流峰值速度曲線與超聲多普勒測速結果十分接近。本實驗處理10只大鼠的造影圖像序列，結果使用超聲粒子成像測速方法得到的峰值速度、平均速度分別較超聲多普勒測速結果低5%~10%和2%~8%。
圖5
超聲粒子成像測速方法與超聲多普勒測速方法的對比
(a)大鼠LCA多普勒測速圖；(b)超聲粒子成像測速法測得最大速度(紅線)與相應區域超聲多普勒測速(綠線)對比圖
Figure5. Ultrasonic particle image velocimetry in comparison with ultrasound Doppler velocimetry(a) blood flow velocity of rat LCA measured by ultrasound Doppler; (b) maximum velocity (red line) measured by ultrasonic particle image velocimetry compared with ultrasound Doppler velocimetry (green line)
通過超聲粒子成像測速技術我們可以獲取任意區域內的血流速度。與超聲多普勒測速相比，超聲粒子成像測速法提供了二維全流場速度信息，能顯示不同位置的血流速度大小和方向，測速效果更為精細全面，能更好地分析流場速度分布情況。通過結合超聲粒子成像測速技術和多項式擬合流速的方法，可計算出全流場剪切率的分布。采用四次多項式擬合血流速度曲線有較好的效果(相關指數R=0.92)，根據擬合曲線計算血管剪切力分布，如圖 6所示。血管中心區域剪切力最小，越靠近血管壁兩邊剪切力越大。
圖6
利用超聲粒子成像測速技術和多項式擬合的方法計算剪切力
(a)多項式擬合的超聲粒子成像測得的血流速度(藍線)；(b)擬合后剪切率分布彩圖
Figure6. Shear calculation by ultrasonic particle image velocimetry combined with polynomial fitting(a) blood flow velocity by ultrasonic particle image velocimetry with polynomial fitting (blue line); (b) shear distribution after polynomial fitting
與利用超聲多普勒測速求剪切力的傳統方法相比，這種新的血管剪切力測量方法使得血管剪切力的計算不單單依靠血管某處速度而是與全流場速度相關[9]，能更有效地評估二維流場剪切力分布。從10只大鼠中選取5只大鼠，使用超聲多普勒測量頸動脈不同位置的血流速度，根據平均血流速度分為低速區、中低速區、中速區、中高速區和高速區，應用本文的方法計算該位置的剪切力，統計三個周期內(約120幀圖)剪切率的平均值，如表 1所示。由于剪切力與剪切率和血液黏度有關，在血液黏度一定的情況下，剪切率變化可反映剪切力的變化。實驗結果表明，大鼠頸動脈血管在周期內受到的平均剪切力與該位置的平均血流速度成正相關。
表1
剪切率與大鼠動脈血流速度的關系
Table1.
Relationship of rat arterial blood flow velocity to shear rate
3 討論與總結
動脈管腔局部狹窄不僅使動脈幾何形狀及力學性質發生顯著改變，更嚴重的是引起動脈中的血液流動異常，誘導動脈組織增生，促使動脈進一步狹窄和粥樣硬化，若斑塊破裂誘發血栓形成，可導致嚴重的后果[10-11]。超聲粒子圖像測速法能動態監測正常血管及狹窄血管的全流場血流動力學變化情況[12-14]，通過提供二維的血流速度和剪切力分布反映血管的狹窄程度和斑塊生長風險，可為臨床提供新的治療依據。本實驗采用頻率高達40 MHz的超聲系統獲取大鼠微小血管和快速心跳頻率下的造影血流[15]；在圖像后處理方面，在圖像互相關算法基礎上結合了迭代算法、亞像素邊緣檢測方法、濾波和插值的偽矢量消除算法；獲得血流二維流場速度分布后，通過多項式擬合的方法計算血流剪切率分布。根據相關研究[16]，采用血管壁附近的血流速度梯度計算剪切力的結果較理論值偏大，而本文通過多項式擬合的方法，模擬流體拋物線的速度分布，使得剪切力計算結果更接近理論值。實驗結果表明利用超聲粒子成像測速可以準確測量復雜血流流場的血流速度和剪切率分布。血管剪切力變化對于心血管疾病的診治具有重要的臨床意義，血管狹窄與斑塊的破裂與長期受到剪切力作用密切相關，但關于血管斑塊生長與剪切力的相互作用機制目前尚未明確，超聲粒子成像測速技術將為進一步測量血管剪切力變化提供一種新的影像方法。
引言
近年來心血管疾病的發病率和死亡率逐年提升，占人類非傳染性疾病總死因的第一位，嚴重威脅著人類健康[1-2]。動脈血管內血流動力學的改變會影響血管壁內皮細胞的形態和功能[2]。究其根本，是由于血液的流動會對血管壁產生包括切應力(剪切力)、周應力和壓應力等在內的各種作用力，從而對血管內皮細胞進行不斷的刺激。其中切應力被認為是與動脈粥樣硬化發生及斑塊破裂最為相關的力學因素[3]。測量血流剪切力要首先測量血管內血流流速分布，并通過剪切力與血流流速梯度之間的正比關系而求出。目前，測量血管血流速度常見的影像學方法是超聲多普勒和磁共振成像(magnetic resonance imaging,MRI)。超聲多普勒只能提供一維的速度信息。為了測量剪切力，需假設血管內流速是標準的層流分布且流速具備拋物線狀分布特征。而這種流速特征只有在形狀規則的血管內才能成立。對于復雜形狀(彎曲、分叉)的血管內的流速分布就會有較大的誤差。而且超聲多普勒技術對聲束與速度向量之間的夾角依賴性大，亦可能產生一定的誤差。MRI具有較高的空間分辨率,但時間分辨率較低，且檢查費用昂貴[4-5]。
超聲粒子成像測速是一種基于粒子成像測速法提出的實時、動態、無創、高分辨率的超聲測速方法。它可以得到精細的二維流場速度矢量、速度梯度分布，進而計算流體剪切力分布[5]。相關研究指出[6]，應用超聲粒子成像測速方法能夠精確計算血管仿體中流體流場的速度分布以及大鼠動脈和靜脈的血流速度。
本文應用超聲粒子成像測速技術獲取大鼠左側頸總動脈(left carotid artery,LCA)血流速度矢量分布，并結合多項式擬合方法計算血管剪切力分布。采集10只大鼠的周期性動脈血流信息，并將超聲粒子成像測速結果與多普勒測速結果進行對比。統計動脈血管不同位置的血流速度和剪切力周期性變化，對于研究血管壁內膜增生和斑塊生長、破裂的機制有重要的意義。
1 材料與方法
超聲粒子圖像測速技術是利用造影微泡對超聲波的強烈反射，并通過高幀頻超聲系統對造影區域采集圖像序列，然后對連續的兩幅圖像進行圖像相關處理，計算圖像中微泡的位移，進而獲得流場的速度和剪切力分布信息。如圖 1所示，超聲粒子圖像測速系統主要由動物超聲系統、圖像序列存儲器、圖像處理算法組成。
圖1
超聲粒子成像測速系統
Figure1.
Ultrasonic particle imaging velocimetry system
1.1 超聲系統
實驗采用由加拿大VisualSonics公司生產的Vevo2100超聲實時分子影像系統。與傳統的超聲系統相比，Vevo2100具備超高的空間分辨率(高達30 μm)和時間分辨率(大于1 000幀/s)。該系統采用了超高頻率的電子線陣探頭(探頭型號為MS400，頻率高達40 MHz)，適用于觀察大鼠微小血管和快速心跳頻率下的造影血流。
1.2 實驗動物及圖像數據采集
10只健康雄性SD大鼠,體重(340±10)g。大鼠被麻醉后，去除左側頸部毛發并固定在試驗臺，使用Vevo2100超聲實時分子影像系統(探頭頻率為40 MHz，幀頻為200 Hz)檢測大鼠左側頸總動脈。超聲檢查始終由同一醫師操作，檢查時超高頻探頭使用支架固定，測量多普勒血流速度。使用5 mL生理鹽水配制超聲造影劑(六氟化碳微泡，聲諾維，SonoVue)并充分搖勻。用1 mL注射器抽取0.2 mL超聲造影劑，經大鼠右側股靜脈注入，再用0.1 mL生理鹽水沖管。同時使用Vevo2100超聲實時分子影像系統對大鼠在造影劑注射時、注射5 min后、注射10 min后分別采集1 000幀圖像，如圖 2所示。
圖2
大鼠LCA的B超圖像
超聲造影劑
Figure2. B-mode images of the rat LCAUCAs: ultrasound contrast agents
1.3 超聲粒子成像測速算法
本文所應用的超聲粒子成像測速法把圖像互相關算法與多次迭代算法、亞像素邊緣檢測方法相結合，并通過濾波和插值的方法消除偽向量，從而改善傳統超聲粒子成像測速算法的測量精度和效率[7]。首先輸入連續的兩幅圖像，設置迭代次數K，再根據圖像互相關算法和亞像素峰值擬合方法計算位移矢量，經插值濾波后得到第一次迭代的結果。若迭代沒有完成，根據前一次迭代的結果進行詢問窗口變形，再重復上述迭代過程。迭代完成后減小詢問窗口大小，應用圖像互相關算法計算位移矢量，然后利用多項式擬合的方法求剪切力，算法流程如圖 3所示。
圖3
算法流程圖
Figure3.
Flow chart of the algorithm
1.3.1 圖像互相關算法
為了獲得流場速度分布，利用圖像互相關算法處理記錄的連續兩幅超聲粒子圖像。首先選擇合適的詢問窗口尺寸，然后超聲粒子圖像被分割成若干個小的詢問窗口，利用互相關函數來衡量兩個詢問窗口之間的灰度相似程度。根據相關文獻的方法[8]，互相關平面內峰值位置直接反映了詢問窗口的位移，再根據連續兩幅圖像的時間差和圖像放大倍數，可獲得該詢問窗口內的流體速度向量。對每一個詢問窗口重復以上步驟就可以獲得二維全流場速度矢量圖。以下是圖像互相關算法公式(1)：
| ${{R}_{ab}}(i,j)=\frac{\sum\limits_{x=0}^{M}{\sum\limits_{y=0}^{N}{[a(i,j)-][b(x-i,y-j)-]}}}{{{\sigma }_{a}}{{\sigma }_{b}}}~,$ |
其中M和N是詢問窗口的維數。a(i,j)和b(i,j)表示灰度值，和表示a和b的平均灰度值，σa和σb是標準差。但是圖像互相關方法在測量亞像素的位移以及粒子團在高速度梯度下的位移時誤差較大，因此本文使用亞像素邊緣檢測方法和多次迭代的方法彌補圖像互相關算法的不足。
1.3.2 亞像素邊緣檢測方法
圖像的灰度是以離散函數的像素單元組成，而位移向量也是以像素為單位。亞像素方法可以評估亞像素的位移，從而提高位移的測量精度。亞像素高斯峰值擬合方法是亞像素分析中常用的方法，以下是亞像素高斯峰值擬合的公式(2)：
| $\left\{ \begin{align} & \Delta x=i+\frac{ln(R(i-1,j))-ln(R(i+1,j))}{2ln(R(i-1,j))-4ln(R(i,j))+2ln(R(i+1,j))} \\ & \Delta y=j+\frac{ln(R(i,j-1))-ln(R(i,j+1))}{2ln(R(i,j-1))-4ln(R(i,j))+2ln(R(i,j+1))} \\ \end{align} \right.,$ |
其中R是互相關算法函數，(i,j)是互相關函數R取得最大值對應的坐標。
1.3.3 多次迭代算法
詢問窗口的粒子團在高速度梯度下流動會產生幾何形變，需要通過詢問窗口變形的方法解決。運行K次迭代算法，第K次迭代時詢問窗口的變形根據第K-1次迭代計算得的位移場Uk-1來獲得。而詢問窗口的變形是通過利用圖像像素點的位移來實現，對于互相關窗口的兩個圖像灰度函數a(i,j)與b(i,j)，利用算法重建后:
| $a\left( r \right)=ar-\frac{{{U}_{K-1}}}{2}-\frac{U{{\prime }_{k-1}}r}{2},$ |
| $b\left( r \right)=br+\frac{{{U}_{K-1}}}{2}-\frac{U{{\prime }_{k-1}}r}{2},$ |
其中r=(i,j),位移梯度U可通過中心有限差分獲得。對經過移位后的圖像再次利用互相關算法和亞像素高斯峰值擬合法計算函數R取得最大值所對應的坐標Rmax，而迭代后的位移UK為：
| ${{U}_{K}}={{U}_{K-1}}+{{R}_{max}}$ |
1.4 多項式擬合法計算血管剪切力
血液在血管腔內流動時，血液與血管壁接觸面上因相對運動產生力F，使血液產生速度變化du。由于液體的粘性將此力層層傳遞，各層液體也相應運動，形成速度梯度du/dr，稱剪切率。設血液與血管壁的接觸面積為A，F/A稱為剪切應力。剪切率與剪切應力具有如下關系：
| $\left( F/A \right)=\eta (du/dr)$ |
其中dr是層厚，η為液體黏度。血管壁受到的剪切力與血流的剪切率成正比，因此利用超聲粒子成像測速計算血流的速度梯度就可評估血管受到的剪切力大小。
血管剪切力是根據血管血流流速具備拋物線狀分布特征來計算的，先應用多項式擬合法擬合拋物線狀分布的血流流速分布，再根據擬合后結果計算剪切力。使用最小二乘法擬合多項式，假設給定數據點(xi,yi)(i=0,1,…,m)，Φ為所有次數不超過n(n≤m)的多項式構成的函數類，求，使得
| $I=\sum\limits_{i=0}^{m}{{{[{{p}_{n}}({{x}_{i}})-{{y}_{i}}]}^{2}}}=\sum\limits_{i=0}^{m}{\left( \sum\limits_{i=0}^{N}{{{a}_{k}}{{x}_{i}}^{k}-{{y}_{i}}^{2}} \right)}=min,$ |
再根據血流流速分布情況確定擬合多項式次數，滿足式(7)的pn(x)稱為最小二乘擬合多項式。
2 結果
本算法基于MATLAB R2009a調試運行，使用24×10至48×16像素大小的詢問窗，0.5~0.75的窗口重疊率，每只大鼠選取其血管微泡濃度適當的圖像序列(約200幀)進行處理。如圖 4所示,選取其中一只大鼠的造影圖像序列，圖 4(a)中黃色虛線區域是大鼠左側頸動脈,在圖像中截取感興趣區域(region of interest,ROI)，選擇像素大小為24×10的詢問窗和50%的窗口重疊率，計算該區域的血流速度分布。圖 4(b)中紅色箭頭的方向和長度分別表示大鼠頸動脈的血流速度方向和大小，并疊加在其對應的血管血流影像區域。如圖 4(c)所示，在血管中心區域血流速度較大，越靠近管壁兩邊血流速度越小。
圖4
應用超聲粒子成像測速技術的大鼠LCA二維血流速度圖
(a)管腔內注射UCAs的大鼠LCA的B超圖像，紅色方框是LCA所在的ROI；(b)ROI的速度矢量分布圖并疊加在其B超圖像上；(c)ROI血流速度梯度分布圖
Figure4. Two-dimensional flow velocity map of rat LCA calculated by ultrasonic particle image velocimetry(a) the B-mode images of the rat LCA with UCAs,and the red box marked the ROI of LCA; (b) velocity vectors map of ROI in the B-mode images; (c) flow velocity gradient map of ROI
為了驗證超聲粒子成像測速方法的可靠性，對連續三個心動周期(約200 ms)超聲粒子成像測速的最大速度取平均值，獲得一個周期的超聲粒子成像測速曲線，并與超聲多普勒測速對比。如圖 5所示，應用本文方法處理其中一只大鼠的數據并與多普勒測速結果對比，超聲粒子成像測速血流峰值速度曲線與超聲多普勒測速結果十分接近。本實驗處理10只大鼠的造影圖像序列，結果使用超聲粒子成像測速方法得到的峰值速度、平均速度分別較超聲多普勒測速結果低5%~10%和2%~8%。
圖5
超聲粒子成像測速方法與超聲多普勒測速方法的對比
(a)大鼠LCA多普勒測速圖；(b)超聲粒子成像測速法測得最大速度(紅線)與相應區域超聲多普勒測速(綠線)對比圖
Figure5. Ultrasonic particle image velocimetry in comparison with ultrasound Doppler velocimetry(a) blood flow velocity of rat LCA measured by ultrasound Doppler; (b) maximum velocity (red line) measured by ultrasonic particle image velocimetry compared with ultrasound Doppler velocimetry (green line)
通過超聲粒子成像測速技術我們可以獲取任意區域內的血流速度。與超聲多普勒測速相比，超聲粒子成像測速法提供了二維全流場速度信息，能顯示不同位置的血流速度大小和方向，測速效果更為精細全面，能更好地分析流場速度分布情況。通過結合超聲粒子成像測速技術和多項式擬合流速的方法，可計算出全流場剪切率的分布。采用四次多項式擬合血流速度曲線有較好的效果(相關指數R=0.92)，根據擬合曲線計算血管剪切力分布，如圖 6所示。血管中心區域剪切力最小，越靠近血管壁兩邊剪切力越大。
圖6
利用超聲粒子成像測速技術和多項式擬合的方法計算剪切力
(a)多項式擬合的超聲粒子成像測得的血流速度(藍線)；(b)擬合后剪切率分布彩圖
Figure6. Shear calculation by ultrasonic particle image velocimetry combined with polynomial fitting(a) blood flow velocity by ultrasonic particle image velocimetry with polynomial fitting (blue line); (b) shear distribution after polynomial fitting
與利用超聲多普勒測速求剪切力的傳統方法相比，這種新的血管剪切力測量方法使得血管剪切力的計算不單單依靠血管某處速度而是與全流場速度相關[9]，能更有效地評估二維流場剪切力分布。從10只大鼠中選取5只大鼠，使用超聲多普勒測量頸動脈不同位置的血流速度，根據平均血流速度分為低速區、中低速區、中速區、中高速區和高速區，應用本文的方法計算該位置的剪切力，統計三個周期內(約120幀圖)剪切率的平均值，如表 1所示。由于剪切力與剪切率和血液黏度有關，在血液黏度一定的情況下，剪切率變化可反映剪切力的變化。實驗結果表明，大鼠頸動脈血管在周期內受到的平均剪切力與該位置的平均血流速度成正相關。
表1
剪切率與大鼠動脈血流速度的關系
Table1.
Relationship of rat arterial blood flow velocity to shear rate
3 討論與總結
動脈管腔局部狹窄不僅使動脈幾何形狀及力學性質發生顯著改變，更嚴重的是引起動脈中的血液流動異常，誘導動脈組織增生，促使動脈進一步狹窄和粥樣硬化，若斑塊破裂誘發血栓形成，可導致嚴重的后果[10-11]。超聲粒子圖像測速法能動態監測正常血管及狹窄血管的全流場血流動力學變化情況[12-14]，通過提供二維的血流速度和剪切力分布反映血管的狹窄程度和斑塊生長風險，可為臨床提供新的治療依據。本實驗采用頻率高達40 MHz的超聲系統獲取大鼠微小血管和快速心跳頻率下的造影血流[15]；在圖像后處理方面，在圖像互相關算法基礎上結合了迭代算法、亞像素邊緣檢測方法、濾波和插值的偽矢量消除算法；獲得血流二維流場速度分布后，通過多項式擬合的方法計算血流剪切率分布。根據相關研究[16]，采用血管壁附近的血流速度梯度計算剪切力的結果較理論值偏大，而本文通過多項式擬合的方法，模擬流體拋物線的速度分布，使得剪切力計算結果更接近理論值。實驗結果表明利用超聲粒子成像測速可以準確測量復雜血流流場的血流速度和剪切率分布。血管剪切力變化對于心血管疾病的診治具有重要的臨床意義，血管狹窄與斑塊的破裂與長期受到剪切力作用密切相關，但關于血管斑塊生長與剪切力的相互作用機制目前尚未明確，超聲粒子成像測速技術將為進一步測量血管剪切力變化提供一種新的影像方法。
