基于矢量場預估的二維超聲彈性成像方法

周佳璐，王韜，楊軍

（中國醫學科學院北京協和醫學院生物醫學工程研究所，天津300192）

1 引 言

傳統臨床中對于軟組織中的囊腫，脂肪瘤等存在彈性差異的病變檢測一般采用觸診的檢測方法，即對組織施加低頻壓力，通過觸覺對病變組織進行定性測量與估計。由于該方法簡單易用性，被廣泛應用于皮膚，乳房，前列腺等部位疾病的檢測中。但觸診本身的客觀性與準確性則有很大的局限，特別是對于病變較深區域或彈性差異較小的病變組織很難進行測量［1］。

準靜態超聲彈性成像是近年來發展較快的一種新的超聲成像技術。其通過對組織施加外部壓力得到組織壓縮前后的超聲回波信號，并結合數字信號處理或者數字圖像技術獲取組織內部不同空間場上的應變信息，從而間接反映組織彈性模量等力學屬性的差異［2］。一方面可以增強傳統超聲成像系統的成像精度，一方面也可以獲取傳統方式獲取不到的病變信息。

超聲彈性成像方法的主要流程一般如下：（1）施加外力獲取前后對比信號；（2）求取圖像位移場；（3）對位移場進行差分求取應變場獲取組織應變［3］。最早的超聲彈性成像法主要采用一維成像方法，即只考慮組織受壓后的軸向位移（沿聲波傳播方向），其主要原因在于求解組織位移場時二維成像的計算量較大，且易引入較大的計算誤差，無法滿足算法的實時性要求，臨床應用性不高［4］。然而組織在受外力發生壓縮形變時，其應變特性是二維的，如果僅考慮軸向位移，計算得到的結果并不全面，精度不高。近年來出現的二維成像方法如動態規劃算法（DP）［5］，快速互相關算法（FNCC）［6］，塊匹配預估算法（BMAPS）［7］，光流法（OF）［8］等同時考慮了彈性形變時的橫向位移，對彈性超聲的臨床應用有一定的推動作用。這些算法分為兩類，一類是基于位移預估的后匹配算法，如塊匹配預估算法（BMAPS）；一類是基于快速匹配方式的算法，如快速互相關算法等。基于快速匹配方式的算法其匹配效果較差，精度較低，臨床使用較少。而現有的基于位移估計的算法在泊松比在0.5附近的情況下效果較好，即其橫向位移不可過大，在泊松比小于0.45的情況下很難獲得較好的結果。

本研究提出了一種基于矢量場預估彈性超聲二維成像方法，采用光流法（OF）進行矢量場預估，采用相位相關算法（PS）進行位移估計，最后采用最小二乘法求取應變值并實現二維成像。實驗采用有限元分析法進行數據仿真與結果對比，并通過對不同泊松比的仿真結果對比證明了該算法可實現較低泊松比情況下成像的清晰穩定。

2 方法

首先獲取包絡線特征峰值指紋信息，采用改進的光流跟蹤法對特征指紋進行跟蹤并進行矢量場預估；然后根據預估矢量場采用相位相關算法更改所跟蹤的RF線對位移場進行精確計算；最后通過最小二乘法濾波求取應變場實現準靜態二維彈性超聲成像。從而可以獲得較高精度的位移場并保證了大尺度橫向應變。

2.1 跟蹤特征獲取

在圖像跟蹤中，一般采用角點跟蹤法，以圖像角點作為光流特征點進行光流跟蹤，這種方法一般適用于分辨率清晰，圖像輪廓分明的普通圖像［9］。而超聲獲取到的射頻信號（RF）歸一化二維圖像特征在于紋理復雜，輪廓并不清晰難以獲取相關角點。因此，提出一種基于包絡線峰值的特征指紋信息作為光流特征點的方法進行光流跟蹤。

首先對RF信號進行歸一化處理，得到RFnorm，有：

其中：min（RF）為射頻信號最小值，max（RF）為射頻信號最大值。

然后求取各條RF信號的包絡線C，有包絡線C滿足：

聯立方程解得上包絡線信號Cup與下包絡線信號Cdown，將RF信號進行分區劃塊，取各區塊RF線上下包絡線峰值與谷值作為特征跟蹤點形成特征指紋信息進行光流跟蹤，見圖1，左圖為立體圖，右圖為俯視圖，橫坐標表示第1～5條RF信號，縱坐標表示信號個數，各點顏色表示包絡信號的波峰波谷位置。

2.2 改進的光流跟蹤

由于RF信號具有紋理特性，如果直接采用傳統光流法，其跟蹤結果置信度不高，因此提出了一種將光流跟蹤與指紋判別檢測算法相結合的方式對光流跟蹤結果進行判斷以去除無效跟蹤信息。

首先，采用光流法對特征點劃分區域進行跟蹤，設各線RF信號包絡線受壓后峰值A變化較小且相對位置特性不變。

設第x條RF線，第y點散射點（x，y）在t時刻的像素值（RF信號幅值）為 A（x，y，t），設壓縮后該點信號幅值表示為 A（x+dx，y+dy，t+dt），o（n）為2階無窮小，則有：

假設各特征點局部區域內包絡信號值恒定不變，即有在各特征點區域內，使得下式O最小即可獲得該特征點區域速度向量

圖1 采用包絡線峰值作為光流跟蹤特征點Fig 1 The envelope peak value is used as the feature points of optical flow tracking

公式中，Ax1，Ay1，At1分別表示跟蹤點（xt，yt）幅值在x，y，t軸的變化率。

然后采用之前作出的特征點矩陣作為指紋特征，對跟蹤結果進行匹配判斷，去除不符合矩陣排布特征的跟蹤結果或置信度1/O較低的結果，其平均值為該區塊的整體速度向量，則矢量跟蹤場結果為：

其中，i11，i12，imn為各區塊中有效值的個數。

2.3 相位相關與應變計算

根據光流跟蹤預估偏移結果（u，v）進行偏移。由于u=dx/dt，壓縮前后只有兩幀數據，則有dt=1，可得u=dx即光流跟蹤時橫向像素偏移量；同理v=dv/dt，則有v=dy即光流跟蹤時縱向像素偏移量。在光流跟蹤處理中，橫向像素數等同于RF線總數，縱向像素數目則是超聲RF信號總采樣點個數，設采樣頻率為fs，則有縱向像素偏移量轉換為時間軸偏移量如下dt=dy/fs=v/s：。設原信號為RF，壓縮后信號為RF。設壓縮前某段信號為RFi（t1，t2），代表第i條RF信號t1到t2之間的一段信號，則壓縮后，該段信號對應的進行相位相關精確匹配的壓縮后的信號段應為 RF（i+u）（t1+v）/f2，t2+v/fs）。

圖2 基于特征指紋的光流矢量跟蹤場示意圖Fig 2 Sketch of optical flow vector tracking field based on feature fingerprint

通過希爾伯特變換對壓縮前后RF信號進行處理，獲取含有相位信息的解析信號，有：

然后，利用跟蹤矢量場作為相關計算的先驗信息，設聲速為s，取預估矢量場橫向位移m=Mu（n，t），縱向時延 tP=Mv（n+m，t）/s，則有壓縮前后 RF解析信號對應的預估偏移信號為：

超聲頻率ω0已知，則通過相關函數相移角度計算兩端信號之間的時延：

則有最終計算得到各時間點t橫向位移為△x=Mu（n，t），縱向位移為△y=s△t。

采用最小二乘法進行應變計算［10］，采用分段線性擬合算法對位移場求差分進行應變擬合即可求得應變場。

3 實驗

3.1 仿真實驗設計

采用有限元仿真軟件COMSOL對模型應變情況進行仿真，采用FIELD II聲場仿真軟件對超聲數據進行仿真處理［11］。

首先用COMSOL建立一個邊長為80 mm的正方形二維模型，其內部中心包含一個半徑為10 mm的圓。設置正方形模型材料與人類組織彈性纖維參數一致：楊氏模量為2MPa，泊松比為0.45。設置圓形材料為病變組織，楊氏模量設為正方形模型材料的5倍：楊氏模量為10MPa，泊松比為0.43［12］。設計受壓模型即正方形底部為有限元分析邊界固定條件，正方形上部為剛性連接，受向下的均衡力為0.1 N/m2。

有限元分析采用Delaunay方法對模型進行自由三角剖分，共產生40000個自由度擴散點，其仿真結果見圖3：

圖3 （a）.有限元分析模型；（b）.模型受壓后Y軸位移場；（c）.模型受壓后X軸位移場Fig 3 （a）.Finite element analysis model；（b）.Y axial displacement field after compression；（c）.X axial displacement field after compression

采用FIELD II聲場仿真軟件設計超聲探頭，其中心頻率為3.5 MHz，采樣率為100 MHz，線陣超聲探頭共計192個陣元，每組發射64個采樣有效陣元，陣元寬度為聲波波長，陣元間隙為波長1/20，（當陣元間隙與陣元寬度之和大于1/2波長時即可避免陣元間相互作用）。陣元長度設置為5 mm，并將每個物理陣元劃分為10個計算陣元。設計聲場掃描模型內部分布40 000個隨機散射點，并根據模型有限元分析結果建立聲場仿真模型見圖4。

圖4 （a）.壓縮前聲場模型；（b）.壓縮后聲場模型Fig 4 （a）.sound field modelbeforecompression；（b）.sound field model after compression

3.2 仿真實驗結果與對比

將矢量場預估算法計算得到的縱向應變結果與有限元分析結果進行對比，比較其偏差值作為算法評價依據，見圖5。

對于傳統一維彈性超聲成像算法，BMAPS塊匹配預估算法，矢量場預估二維成像算法在0.4～0.48不同泊松比情況下進行對比，結果見表1。

由結果可看出，傳統一維彈性超聲算法在較低泊松比情況下其計算誤差成指數增長，經過數據觀察，其只能在中心RF線上保證較好的計算結果，而本算法在較低泊松比條件下的計算誤差也低于BMAPS塊匹配預估算法。實驗證明，本算法是一種清晰度更高，魯棒性更好的超聲彈性成像方法。

表1 不同泊松比情況下各算法的偏差性能對比Table 1 Comparison of the deviation performance of each algorithm under Different Poisson's ratio

圖5 （a）應變計算結果（b）有限元分析計算結果Fig 5 （a）strain calculation results（b）finite element analysis and calculation results

4 總結

本研究旨在提出一種高清晰度高魯棒性的超聲彈性成像方法，采用光流跟蹤法對包絡線峰值指紋特征進行跟蹤，以跟蹤矢量場結果作為先驗信息對超聲RF信號進行相位相關計算，一方面實現橫向位移的準確估計，一方面也避免了相位相關計算時的相位混疊問題。與傳統一維相位相關法相比，本算法在計算相關函數之前將跟蹤矢量場位移作為預估量加入計算，避免在存在橫向位移時依舊同一條RF線上進行信號匹配，與近年來的二維相位相關法相比，采用光流法預估矢量場大大提高了橫向位移預估的精確性，可保證在軟組織泊松比較小的情況下依舊可以進行準確追蹤。實驗結果最終也證明了該算法的準確有效性，對于超聲彈性成像技術的系統研發具有一定的推動作用。