任雅姿， 李　穎 ， 劉　歡

(河北工業(yè)大學(xué) 電氣工程學(xué)院，天津　300130)

?

基于Hankel矩陣的奇異值分解法對腦血容積計算的研究

任雅姿， 李穎 ， 劉歡

(河北工業(yè)大學(xué) 電氣工程學(xué)院，天津300130)

針對動態(tài)對比度增強磁共振灌注成像中腦血容積的計算，提出基于Hankel矩陣的奇異值分解(Singular Value Decomposition，SVD)算法。在奇異值數(shù)目的確定上采用差分譜量級差的研究方法，對算法進行理論推導(dǎo)與仿真模擬，得到較為理想的濾波效果。由于成像過程存在測量噪聲的干擾，分析了信噪比和示蹤劑延遲對算法的影響。仿真結(jié)果表明，信噪比越低(SNR=5 dB)，算法處理效果越明顯；信噪比增高(SNR=100 dB)，估計值偏差減小，結(jié)果越為準確。且該算法不受示蹤劑延遲的影響。與傳統(tǒng)奇異值分解算法相比，采用基于Hankel矩陣的奇異值算法可以更為準確地估計腦血容積。

動態(tài)對比度增強磁共振灌注成像；腦血容積；奇異值分解；Hankel矩陣

動態(tài)對比度增強磁共振灌注成像技術(shù)在估計腦缺血和腦腫瘤微循環(huán)血液動力學(xué)的變化方面有直觀準確性，可以提供組織動態(tài)增強信息，且能無創(chuàng)地反映腫瘤組織局部血流灌注和毛細血管通透性等血流動力學(xué)方面的信息，如腦血容積(Cerebral Blood Volume，CBV)、腦血流量(Cerebral Blood Flow，CBF)、平均通過時間(Mean Transit Time，MTT)等[1]。

由于腦血容積求取過程以及所需相關(guān)信號的特殊性，對所需信號的濾波處理方法也受到限制，目前主要是傅里葉變換法和奇異值分解法[2]。研究表明，傅里葉變換法對噪聲較敏感，在信噪比低時腦血容積有一定程度的低估；而奇異值分解法受閾值限制，處理過程較為復(fù)雜[3-5]。為使處理過程簡便，結(jié)果客觀，本文提出一種不受閾值影響的新算法——基于Hankel矩陣的奇異值分解法。在奇異值數(shù)目的確定上，采用基于差分譜的量級差的改進方法，結(jié)果表明，該方法能夠準確選取與最優(yōu)降噪效果相對應(yīng)的奇異值數(shù)目，在對含噪信號的降噪處理上顯示了良好的性能。進一步，應(yīng)用改進的奇異值分解法進行了信噪比和示蹤劑延遲等影響因素的分析。

1　腦血容積計算的理論推導(dǎo)

1.1基本理論

臨床上，動態(tài)對比度增強的核磁共振灌注成像診斷疾病所用的血液動力學(xué)參數(shù)都是基于示蹤劑濃度的。示蹤劑濃度與腦血流量(CBF)的關(guān)系式[6]如下：

CBF·Ca(t)?R(t)

(1)

式中:Ct(t)為示蹤劑濃度，Ca(t)為動脈輸入函數(shù)(Arterial Input Function，AIF)，CBF為腦血流量，R(t)為駐留函數(shù)，表示示蹤劑到達組織后隨時間衰減的幾率，為衰減型函數(shù)，且當(dāng)t=0時，R=1。

感興趣區(qū)的腦血容積可用示蹤劑濃度的面積積分與動脈輸入函數(shù)的面積積分之比求出。

(2)

式中:CBV為腦血容積。

在臨床應(yīng)用中，信號強度與示蹤劑濃度關(guān)系函數(shù)：

(3)

式中:S0為示蹤劑未到達組織之前測得的信號強度，TE為回波時間，k為擬合常數(shù)。k值的確定包含了諸多的因素，它與腦組織密度和血細胞比容有關(guān)，在不同的動脈、靜脈和組織中有不同的k值，其大小一般通過實驗的方法確定。

將式(3)變換，可得：

S(t)=S0e-(kCt(t)TE)

(4)

在臨床應(yīng)用中，信號的采集要受到噪聲的干擾。在MR信號中，萊斯噪聲是主要的噪聲，當(dāng)信噪比(SNR)大于2時，萊斯噪聲與高斯白噪聲相似[7]。因此，本實驗中采用高斯白噪聲近似代替萊斯噪聲。這樣，實際測量的信號可以近似表示為：

Sn(t)=S0e-(kCt(t)TE)+N(t)

(5)

將式(5)變換，得到含噪的示蹤劑濃度函數(shù)：

(6)

上述理論得對含噪的示蹤劑濃度進行濾波處理即可得到精確腦血容積值。

1.2利用基于Hankel矩陣的奇異值分解法求腦血容積的基本原理

1.2.1Hankel矩陣與SVD分解

將式(6)表示為：

Cn=[Cn(1),Cn(2),…,Cn(k),…,Cn(N)]

(7)

式中:k=1，2，…，N，其中N為數(shù)據(jù)長度。用上述數(shù)據(jù)構(gòu)造p×q階Hankel矩陣，有：

(8)

式中:p+q-1=N,p≥q,H(i,j)=Cn(i+j-1)。

將構(gòu)造的Hankel矩陣H進行奇異值分解可得到：

H=U∑VT

(9)

式中:U為p×p維正交左奇異陣；V為q×q維正交右奇異陣；∑為p×q維矩陣。主對角元素為

λi(i=1,2,…k),k=min(p,q)

即：

∑=diag(λ1,λ2,…,λk)

(10)

式中:λ1,λ2,…,λk為矩陣H的奇異值，且

λ1≥λ2≥…≥λk≥0

矩陣H是由噪聲污染的信號構(gòu)成Hankel矩陣，其中的有用信號和噪聲信號具有互不相關(guān)性，所以可以表示為互不相關(guān)的兩個子空間，即無噪聲污染的信號子空間Hs和噪聲子空間Hn之和，則有：

(11)

從式(11)可得，對原始信號的降噪就是從Hankel矩陣H中去除噪聲空間Hn，即已知H，尋找Hs的最佳逼近問題，逼近程度越好，降噪效果越明顯。

傳統(tǒng)的SVD分解原理很簡單，不需要構(gòu)造Hankel陣，分解矩陣直接重構(gòu)就可獲得濾波后的信號。但是，對于重構(gòu)之前的閾值的選擇要求相當(dāng)高，也沒有可行的理論選擇法，主要依靠經(jīng)驗和反復(fù)修正，在模擬中可以實現(xiàn)，但是對于真實信號的處理相當(dāng)不理想。

1.2.2重構(gòu)陣維數(shù)與有效奇異值個數(shù)的確定

根據(jù)奇異值分解和Frobeious范數(shù)的矩陣最佳逼近定理可得：有用信號主要由前r個相對較大的奇異值反映，噪聲信號由后面相對較小的奇異值反映，所以去掉代表噪聲信號的較小奇異值，則從原始信號中去除了噪聲。確定有效奇異值數(shù)目后，進行奇異值分解的逆過程，從重構(gòu)矩陣中提取降噪信號。因此，重構(gòu)陣的維數(shù)和有效奇異值數(shù)目直接影響降噪效果。

大量研究實驗表明，Hankel重構(gòu)矩陣的最佳維數(shù)基本在p=N/2處的一個鄰域內(nèi)產(chǎn)生，其降噪效果比較理想，能滿足要求[8]。所以本文的重構(gòu)陣的結(jié)構(gòu)根據(jù)Cn的長度N來確定，取p=N/2。

對于有效奇異值個數(shù)的確定，采用奇異值的差分譜來具體確定。奇異值的差分譜序列ai=λi-λi+1,(i=1,2,…,q-1)描述了奇異值序列的具體變化情況。差分序列A=(a1,a2,…,aq-1)充分顯示了相鄰兩個奇異值的變化。其差別越大，在整個差分譜中表現(xiàn)的特征也越明顯，它們之間的差值也相對越大。較大的值通常會攜帶非常重要的狀態(tài)信息，其出現(xiàn)是由于有用信號和噪聲信號的不相關(guān)性導(dǎo)致的。所以，根據(jù)差值的不同，選取合適的奇異值數(shù)目可達到去噪目的。本文加噪信號中噪聲的奇異值是一條連續(xù)平滑的幅度很小的無拐點曲線，接近于一條平坦的直線。根據(jù)加噪信號的特點，本文采用一種新的確定奇異值數(shù)目的方法，即差分譜從右到左，選擇第一個相鄰兩個值的差增加一個量級時(相鄰兩個值的商大于10)的奇異值的序列數(shù)為有效奇異值數(shù)目k。

1.2.3Hankel陣重構(gòu)與腦血容積估計

Ct(1)=Hs(1,1),

?

?

Ct(N)=Hs(p,q)

(12)

2　腦血容積計算的仿真實驗

針對基于Hankel矩陣的奇異值分解算法計算腦血容積的準確性問題，進行了以下幾個方面的仿真。

2.1腦血容積理論值的構(gòu)造

腦血容積的仿真模擬中，典型的動脈輸入函數(shù)如下[10]：

(13)

式中:C0=1，t0=10 s，r=3，b=1.5，將60 s作為一個完整的動脈輸入函數(shù)。

為驗證算法的準確性，采用如下三種理想的駐留函數(shù)[10]：

(1) 指數(shù)(Exponential)型駐留函數(shù)

(14)

(2) 三角(Triangle)型駐留函數(shù)

(15)

(3) 盒(Box-shaped)型駐留函數(shù)

(16)

式中:MTT為平均通過時間。仿真中的MTT=10 s，CBF=60 mL/100 mg/min。由式(1)與式(2)可得理論腦血容積CBV=10.508 mL/100 mg。

2.2腦血容積模擬值的構(gòu)造

臨床上可檢測組織信號強度，由此可求腦血容積值。其中理想的信號強度與示蹤劑濃度的關(guān)系見式(4)。仿真中對理想信號強度加噪來模擬實際測量信號，其與示蹤劑濃度關(guān)系見式(5)。其中取原始信號強度S0=100，回波時間TE=0.036 s，k值選擇在信號強度曲線下降原來的40%處。信噪比選擇SNR=10 dB。

2.3基于Hankel陣的SVD算法計算腦血容積

采用指數(shù)型駐留函數(shù)仿真，處理前后的示蹤劑濃度與算法的收斂性見圖1。其中，圖1(a) 表示理想示蹤劑濃度，圖1(b)表示加噪后的示蹤劑濃度，圖1(c)表示降噪處理后的示蹤劑濃度，圖1(d)表示奇異值分布曲線。

由圖1可得當(dāng)奇異值數(shù)目為4時，可達到較理想的濾波效果。其中，經(jīng)計算，理想CBV=10.508 mL/100 mg，加噪后的CBV1=13.336 mL/100 mg，算法降噪處理后的CBV=10.938 mL/100 mg。所以采用差分譜量極差的方法確定奇異值數(shù)目的基于Hankel的SVD算法能夠較為準確地估計腦血容積。

3　影響因素的分析

3.1信噪比對算法的影響

在動態(tài)對比度增強磁共振灌注成像過程中會受到不同程度噪聲的影響。在圖2中，分別顯示信噪比SNR=5 dB時，算法降噪處理前后的三種不同駐留函數(shù)下的腦血容積估計值。其中，圖2(a)表示降噪前的CBV值，圖2(b)表示降噪后的CBV值。圖3為信噪比SNR=10 dB時，算法降噪處理前后的三種不同駐留函數(shù)下的腦血容積估計值。其中，圖3(a)表示降噪前的CBV值，圖3(b)表示降噪后的CBV值。圖4為高信噪比下，算法處理后的CBV估計值，其中圖4(a)表示信噪比SNR=50 dB時的CBV值，圖4(b)表示信噪比SNR=100 dB時的CBV值。

圖1　基于Hankel矩陣的SVD算法處理Cn的效果圖Fig.1　Cn processed by Hankel matrix and SVD algorithm

由圖2和圖3可知，不同信噪比對腦血容積值有一定影響，隨著信噪比的增大，降噪處理前所得CBV值高估現(xiàn)象降低?；贖ankel矩陣的SVD算法處理后，與處理前的值相對比，有明顯改善。進一步，該算法在信噪比越小的情況下，降噪效果越明顯。由圖4可知，在高信噪比下，腦血容積值的估計無偏差，經(jīng)過算法處理后的值較為準確。綜合圖2～圖4，可得，在同一信噪比下，不同駐留函數(shù)下的腦血容積估計具有一致性。

3.2示蹤劑延遲時間對算法的影響

示蹤劑一般會先到達正常動脈端，再到達組織，因此示蹤劑延遲在腦血液動力學(xué)參數(shù)的測量中不可避免。模擬示蹤劑延遲對CBV的影響時的數(shù)學(xué)表達式為

Ct(t-Td)=Ca(t)?[CBF·R(t-Td)]

(17)

式中:Td為示蹤劑延遲時間，即對比劑到達動脈端與到達組織的時間差。

為驗證單一因素對腦血容積估計值的影響，采用指數(shù)型駐留函數(shù)，參考CBV分別為2 mL/100 mg和8 mL/100 mg，SNR=100 dB，模擬示蹤劑延遲時間為0～10 s，結(jié)果見圖5，示蹤劑延遲時間對基于Hankel的SVD算法估計腦血容積無影響。

圖5　不同延遲時間下的CBVFig.5 CBV estimation with different delays

4　結(jié)　論

本研究對基于Hankel矩陣SVD算法估計腦血容積進行了理論推導(dǎo)，仿真了三種不同駐留函數(shù)下的腦血容積值的估計，分析得到高信噪比與示蹤劑延遲時間對算法均無影響。采用差分譜量級差的方法確定奇異值數(shù)目，得到了SVD算法可以在高信噪比下準確估計任意駐留函數(shù)下的腦血容積。克服了傳統(tǒng)奇異值數(shù)目選取不當(dāng)對處理結(jié)果影響的弊端，該算法的穩(wěn)定性和準確性為臨床上通過腦血容積診斷疾病提供了參考依據(jù)。

[1] EDELAN R R, WIELOPOLSKI P, SCHMITT F. Echo-planar MR imaging[J]. Radiology, 1994, 193(2): 600-612.[2] OSTERGAARD L, WEISSKOFF R M, CHESLER D A, et al. High resolution measurement of cerebral blood flow using intravascular tracer bolus passages[J]. Magnetic Resonance Imaging, 1996, 36(5): 715-725．

[3] WU O, OSTERGAARD L, KOROSHETZ W J, et al. Effects of tracer arrival time on flow estimates in MR perfusion-weighted imaging[J]. Magnetic Resonance Imaging, 2003, 50(4): 856-864.

[4] SMITH M R, LU H, TROCHET S, et al. Removing the effect of SVD algorithmic artifacts present in quantitative MR perfusion studies[J]. Magnetic Resonance Imaging, 2004, 51(4): 631-634.

[5] LI Y, MA D, HE R, et al. The application of magnetic resonance perfusion imaging in the estimation of brain function using SVD method[J]. IEEE Transactions on Magnetics, 2012, 48(11): 2850-2853.

[6] 俎棟林. 核磁共振成像學(xué)[M].北京：高等教育出版社, 2004.

[7] UNAL S, ROHIT S. Cerebral blood flow estimation from perfusion-weighted MRI using FT-based MMSE filtering method[J].Magnetic Resonance Imaging, 2008, 26(3): 313-322.

[8] 王建國, 李健, 劉穎源. 一種確定奇異值分解降噪有效秩階次的改進方法[J]. 振動與沖擊,2014, 33(12):176-180.

WANG Jianguo, LI Jian, LIU Yingyuan. An improved method for determining effective order rank of SVD denoising [J].Journal of Vibration and Shock, 2014, 33(12):176-180.

[9] 曾作欽. 基于奇異值分解的信號處理方法及其在機械故障診斷中的應(yīng)用[D].廣州：華南理工大學(xué), 2011.

[10] PERUZZO D, BERTOLDO A, ZANDERIGO F, et al. Automatic selection of arterial input function on dynamic contrast-enhanced MR images[J]. Computer Methods and Programs in Biomedicine, 2011, 104(3):148-157.

Calculation of cerebral blood volume using the SVD method based on the Hankel matrix

REN Yazi, LI Ying, LIU Huan

(School of Electrical Engineering, Hebei University of Technology, Tianjin 300130, China)

The singular value decomposition(SVD) method based on the Hankel matrix was proposed in this paper to calculate cerebral blood volume(CBV) in dynamic contrast-enhanced perfusion magnetic resonance imaging. The method of differential spectral order error was proposed to determine the number of singular values. Through the theoretical derivation and simulation, the ideal filtering effect was reached. Because of the existing noise in the imaging process, the influences such as signal-to-noise ratio(SNR) and tracer delay were analyzed. The simulation results show that this method can estimate cerebral blood volume more effectively at the higher SNR(SNR=100 dB) and is not affected by the tracer delay. While compared with the value of the cerebral blood volume without and with using the method, the superiority of the method is clearly shown at the lower SNR(SNR=5 dB). Compared with the traditional singular value decomposition, the SVD method based on the Hankel matrix can be used to estimate the cerebral blood volume more accurately.

dynamic contrast-enhanced perfusion magnetic resonance imaging; cerebral blood volume; singular value decomposition (SVD); Hankel matrix

河北省自然科學(xué)基金項目資助(E2013202238)

2015-07-17修改稿收到日期：2015-08-26

任雅姿 女，碩士生，1989年生

李穎 女，博士，教授，博士生導(dǎo)師，1973年生

R318

A

10.13465/j.cnki.jvs.2016.16.007