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

任雅姿， 李　穎 ， 劉　歡

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

?

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

任雅姿， 李穎 ， 劉歡

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

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

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

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

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

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

1.1基本理論

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

CBF·Ca(t)?R(t)

(1)

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

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

(2)

式中:CBV為腦血容積。

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

(3)

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

將式(3)變換，可得：

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

(4)

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

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

(5)

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

(6)

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

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ù)長(zhǎng)度。用上述數(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進(jìn)行奇異值分解可得到：

H=U∑VT

(9)

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

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

即：

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

(10)

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

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

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

(11)

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

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

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

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

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

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

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

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

?

?

Ct(N)=Hs(p,q)

(12)

2　腦血容積計(jì)算的仿真實(shí)驗(yàn)

針對(duì)基于Hankel矩陣的奇異值分解算法計(jì)算腦血容積的準(zhǔn)確性問(wèn)題，進(jìn)行了以下幾個(gè)方面的仿真。

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

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

(13)

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

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

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

(14)

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

(15)

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

(16)

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

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

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

2.3基于Hankel陣的SVD算法計(jì)算腦血容積

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

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

3　影響因素的分析

3.1信噪比對(duì)算法的影響

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

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

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

3.2示蹤劑延遲時(shí)間對(duì)算法的影響

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

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

(17)

式中:Td為示蹤劑延遲時(shí)間，即對(duì)比劑到達(dá)動(dòng)脈端與到達(dá)組織的時(shí)間差。

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

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

4　結(jié)　論

本研究對(duì)基于Hankel矩陣SVD算法估計(jì)腦血容積進(jìn)行了理論推導(dǎo)，仿真了三種不同駐留函數(shù)下的腦血容積值的估計(jì)，分析得到高信噪比與示蹤劑延遲時(shí)間對(duì)算法均無(wú)影響。采用差分譜量級(jí)差的方法確定奇異值數(shù)目，得到了SVD算法可以在高信噪比下準(zhǔn)確估計(jì)任意駐留函數(shù)下的腦血容積?？朔藗鹘y(tǒng)奇異值數(shù)目選取不當(dāng)對(duì)處理結(jié)果影響的弊端，該算法的穩(wěn)定性和準(zhǔn)確性為臨床上通過(guò)腦血容積診斷疾病提供了參考依據(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] 王建國(guó), 李健, 劉穎源. 一種確定奇異值分解降噪有效秩階次的改進(jìn)方法[J]. 振動(dòng)與沖擊,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] 曾作欽. 基于奇異值分解的信號(hào)處理方法及其在機(jī)械故障診斷中的應(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é)基金項(xiàng)目資助(E2013202238)

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

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

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

R318

A

10.13465/j.cnki.jvs.2016.16.007