周佳璐，王韜，楊軍

（中國醫(yī)學(xué)科學(xué)院北京協(xié)和醫(yī)學(xué)院生物醫(yī)學(xué)工程研究所，天津300192）

1 引 言

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

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

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

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

2 方法

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

2.1 跟蹤特征獲取

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

首先對RF信號進(jìn)行歸一化處理，得到RFnorm，有：

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

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

聯(lián)立方程解得上包絡(luò)線信號Cup與下包絡(luò)線信號Cdown，將RF信號進(jìn)行分區(qū)劃塊，取各區(qū)塊RF線上下包絡(luò)線峰值與谷值作為特征跟蹤點(diǎn)形成特征指紋信息進(jìn)行光流跟蹤，見圖1，左圖為立體圖，右圖為俯視圖，橫坐標(biāo)表示第1～5條RF信號，縱坐標(biāo)表示信號個數(shù)，各點(diǎn)顏色表示包絡(luò)信號的波峰波谷位置。

2.2 改進(jìn)的光流跟蹤

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

首先，采用光流法對特征點(diǎn)劃分區(qū)域進(jìn)行跟蹤，設(shè)各線RF信號包絡(luò)線受壓后峰值A(chǔ)變化較小且相對位置特性不變。

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

假設(shè)各特征點(diǎn)局部區(qū)域內(nèi)包絡(luò)信號值恒定不變，即有在各特征點(diǎn)區(qū)域內(nèi)，使得下式O最小即可獲得該特征點(diǎn)區(qū)域速度向量

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

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

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

其中，i11，i12，imn為各區(qū)塊中有效值的個數(shù)。

2.3 相位相關(guān)與應(yīng)變計(jì)算

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

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

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

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

超聲頻率ω0已知，則通過相關(guān)函數(shù)相移角度計(jì)算兩端信號之間的時延：

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

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

3 實(shí)驗(yàn)

3.1 仿真實(shí)驗(yàn)設(shè)計(jì)

采用有限元仿真軟件COMSOL對模型應(yīng)變情況進(jìn)行仿真，采用FIELD II聲場仿真軟件對超聲數(shù)據(jù)進(jìn)行仿真處理［11］。

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

有限元分析采用Delaunay方法對模型進(jìn)行自由三角剖分，共產(chǎn)生40000個自由度擴(kuò)散點(diǎn)，其仿真結(jié)果見圖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聲場仿真軟件設(shè)計(jì)超聲探頭，其中心頻率為3.5 MHz，采樣率為100 MHz，線陣超聲探頭共計(jì)192個陣元，每組發(fā)射64個采樣有效陣元，陣元寬度為聲波波長，陣元間隙為波長1/20，（當(dāng)陣元間隙與陣元寬度之和大于1/2波長時即可避免陣元間相互作用）。陣元長度設(shè)置為5 mm，并將每個物理陣元劃分為10個計(jì)算陣元。設(shè)計(jì)聲場掃描模型內(nèi)部分布40 000個隨機(jī)散射點(diǎn)，并根據(jù)模型有限元分析結(jié)果建立聲場仿真模型見圖4。

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

3.2 仿真實(shí)驗(yàn)結(jié)果與對比

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

對于傳統(tǒng)一維彈性超聲成像算法，BMAPS塊匹配預(yù)估算法，矢量場預(yù)估二維成像算法在0.4～0.48不同泊松比情況下進(jìn)行對比，結(jié)果見表1。

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

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

圖5 （a）應(yīng)變計(jì)算結(jié)果（b）有限元分析計(jì)算結(jié)果Fig 5 （a）strain calculation results（b）finite element analysis and calculation results

4 總結(jié)

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