奇異值分解法用于MR灌注成像腦血流量估計的仿真研究

馬棟敏 李 穎* 彭 川 郭延慶 郭 磊 何任杰 饒利蕓5

（河北工業(yè)大學(xué)河北省電磁場與電器可靠性省部共建重點實驗室，天津 300130）2（天津市第四中心醫(yī)院，天津 300140）3（軍械工程學(xué)院，石家莊 050003）4（University of Texas Medical School at Houston，Houston，TX77030，USA）5（The Methodist Hospital Research Institute，Houston，TX 77030，USA）

引言

MR腦灌注成像（perfusion imaging）是將組織毛細血管水平的血流灌注情況，通過磁共振成像方式顯示，以評估局部的組織活力及功能（代謝能力）。它通過采用示蹤劑團注技術(shù)的原始圖像取得組織的信號強度-時間曲線和動脈輸入函數(shù)（arterial input function，AIF），提供血液動態(tài)代謝能力信息，如腦血容量（cerebral blood volume，CBV）、腦血流量（cerebral blood flow，CBF）、平均通過時間（mean transit time，MTT）等參數(shù)［1-2］。可分析、評價組織器官的灌注狀態(tài)，得到病變部位組織的微循環(huán)血流及血管變化情況等信息。

MR腦灌注成像參數(shù)的確定是一個逆問題求解的過程。Miles等提出可以用最大斜率法確定灌注參數(shù)［3］，這種方法原理最簡單，但是要求示蹤劑注射速率在8～10 mL/s，過高的對比劑注射速率和注射量不僅增加注射部位對比劑外滲的危險性，且要求病人有良好的心肺功能，因此這種方法在臨床上并不適用。1996年，?stergaard等提出傅里葉變換和奇異值分解法都可以用來估計灌注參數(shù)［4］。研究發(fā)現(xiàn)傅里葉變換對噪聲敏感，在信噪比低時會造成血流量的低估［5-9］，因此奇異值分解法受到人們更廣泛的關(guān)注。在我國，人們對MR灌注成像的研究主要在臨床應(yīng)用方面，對灌注參數(shù)獲取方法的研究少之又少。文中對奇異值分解法估計腦血流量進行了簡單介紹，對不同閾值選擇方式、信噪比、動脈輸入函數(shù)的延遲和失真等因素的影響進行了仿真實驗，針對奇異值分解法在動脈輸入函數(shù)延遲時低估血流量的情形進行了修正，以使該方法更有效。

1 理論推導(dǎo)

采用示蹤劑團注技術(shù)確定CBV、CBF和MTT主要基于示蹤劑稀釋理論。示蹤劑以團注的形式從受試者靜脈進入組織的過程可表示為［4］

式中，?表示卷積運算，ct（t）表示組織示蹤劑濃度，ca（t）表示動脈輸入函數(shù)，R（t）表示駐留函數(shù)（residue function），即示蹤劑到達組織后隨著時間還留在組織內(nèi)的幾率，R（t=0）=1。

將式（1）離散化并寫成矩陣的形式

將式（3）寫為

式中，C表示最左邊的矩陣，B表示最右邊的矩陣，其余的用A表示。

可以從灌注所得原始圖像中得到ct（t）和ca（t），為了得到 CBF，需對矩陣 A求逆，將 A進行奇異值分解，得

有

則

因為在灌注成像的過程中，測量噪聲和生理噪聲不可避免，為了降低噪聲對 CBF的影響，選定閾值Psvd，將對角陣W 中小于Psvd的值給0，然后將W代入式（7）可得B。

2 仿真實驗

為了驗證奇異值分解法計算腦血流量的準確性，并分析其影響因素，針對不同的信噪比及動脈輸入函數(shù)的延遲與失真進行了幾組仿真實驗。仿真試驗包括以下幾個方面。

2.1 模擬動脈輸入函數(shù)（AIF）

所用動脈輸入函數(shù)為對受試者注入標準劑量的示蹤劑所能得到的經(jīng)過伽馬擬合的較為典型的

式中，r=3，b=5，c0=1，t0=10 s，整個曲線過程為60 s。

當(dāng)腦組織發(fā)生病變時，會出現(xiàn)動脈輸入函數(shù)的延遲或失真，這些情況都會對腦血流量的估計產(chǎn)生影響。在仿真實驗中，通過將動脈輸入函數(shù)循環(huán)平移來仿真AIF延遲；通過改變 r、b的值以改變動脈輸入函數(shù)的形狀來仿真AIF失真。

2.2 模擬組織示蹤劑濃度-時間曲線

為了驗證算法的準確性，假設(shè)已經(jīng)知道駐留函數(shù)R（t），在理想的情況下，常用的駐留函數(shù)有3種描述［4］。

1）指數(shù)型（exponential）駐留函數(shù)一種函數(shù)［10］，其數(shù)學(xué)表達式為

2.3 模擬信號強度-時間曲線

組織的信號強度與示蹤劑濃度之間的關(guān)系為［11-12］

組織示蹤劑濃度-時間曲線由式（1）可得，根據(jù)中心容積定理有：

3）盒型（box）駐留函數(shù)

取s0=100，TE=66 ms，ct（t） 轉(zhuǎn)換成st（t） 的k值是選擇在 CBV=4 mL/100 mg，CBF=60 mL/100 mg/min時，信號曲線下降了原來的40％，ca（t） 轉(zhuǎn)換成sa（t） 的k值是選擇在 CBV=4 mL/100 mg，CBF=60 mL/100 mg/min時，信號曲線下降了原來的 60％［10］。

2.4 采用奇異值分解法作為仿真算法

考慮到生理噪聲和測量噪聲的影響，為了更好的模擬實際的信號強度-時間曲線，對st（t） 和sa（t）加上均值為0，標準差為σ的高斯白噪聲，σ與信噪比 SNR 之間的關(guān)系［13］為

將加入噪聲后的信號強度-時間曲線按式（13）變換為示蹤劑濃度-時間曲線，將所得曲線采用奇異值分解法得到腦血流量估計結(jié)果。

將動脈輸入函數(shù)ca（t）生成的矩陣A采用Matlab7.1中自帶的奇異值分解函數(shù)進行奇異值分解得到矩陣U、V和W。奇異值分解法通過將對角陣W中小于某一閾值的奇異值置零來濾除噪聲，所以閾值選擇是否合適直接影響奇異值分解法對腦血流量的估計結(jié)果。

在信噪比為150和10時，對不同駐留函數(shù)不同的血流量值選擇相同的閾值，在低血流量時，估計結(jié)果很好；但在高血流量時，估計結(jié)果變差。針對以上情況，提出對不同駐留函數(shù)和不同血流量值，選擇不同的閾值，在高血流量時選擇比低血流量時小的閾值。不同的閾值水平分別為矩陣W中最大值的15％和矩陣W中最大值的30％。

2.5 延遲效應(yīng)的修正

當(dāng)腦組織發(fā)生病變時，會出現(xiàn)動脈輸入函數(shù)的延遲或失真，這些情況都會對腦血流量的估計產(chǎn)生影響。AIF延遲會造成高血流量的低估，所以提出對延遲效應(yīng)進行修正。

修正的基本思想是：首先對示蹤劑濃度-時間曲線峰值前的部分進行最小二乘擬合以得到延遲時間，然后根據(jù)此延遲時間將組織示蹤劑濃度-時間曲線進行時移，最后利用時移后的曲線采用奇異值分解法得到腦血流量的估計結(jié)果。

因為噪聲的影響，每次計算結(jié)果會有所不同，故在仿真實驗中，每個血流量值計算200次，取均值為最終結(jié)果，且對每個血流量值200次估計結(jié)果進行了標準差的計算。標準差（sd）為

式中，xi為每次血流量估計值，ˉx為均值，n為計算次數(shù)200。

3 結(jié)果

3.1 不同信噪比和閾值水平的仿真結(jié)果

圖1為SNR=150，分別采用3種駐留函數(shù)，閾值分別為矩陣W中最大值的15％和30％時的結(jié)果。

由圖可見，在低血流量時，閾值對腦血流量的估計結(jié)果影響不明顯，估計結(jié)果很好；但在高血流量時，閾值對腦血流量的估計結(jié)果有較明顯的影響，估計結(jié)果變差。這是因為在血流量值大時，其示蹤劑濃度-時間曲線的峰值要大于血流量值小時的情況。顯然，在不同血流量值時選擇相同的閾值，會造成高血流量時有用信息的濾除或低血流量時噪聲不能完全濾除。此外，由圖1還可以看出，駐留函數(shù)為盒型駐留函數(shù)時，在高血流量時會有高估的情形出現(xiàn)。這是因為盒型駐留函數(shù)過于理想，不能很好地描述示蹤劑到達組織后隨著時間還留在組織內(nèi)的幾率。

圖1 SNR=150，不同閾值水平的腦血流值估計結(jié)果。（a）閾值為15％；（b）閾值為30％Fig.1 Estimated CBF values with SNR=150.（a）Psvd=15％；（b）Psvd=30％

圖2為SNR=10，分別采用3種駐留函數(shù)，閾值分別為矩陣W中最大值的15％和30％時的結(jié)果。由結(jié)果可以看出，信噪比較高時，腦血流量估計較準確，且每次計算結(jié)果比較集中，此時該方法估計腦血流量結(jié)果比較穩(wěn)定。當(dāng)SNR較大時，閾值對腦血流量較大時的估計結(jié)果有較明顯的影響；當(dāng)SNR較小時，為了估計結(jié)果比較準確，需要較大的閾值來濾除噪聲的影響。

3.2 不同閾值的影響

圖2 SNR=10，不同閾值水平的腦血流值估計結(jié)果。（a）閾值為15％；（b）閾值為30％Fig.2 Estimated CBF values with SNR=10.（a）Psvd=15％；（b）Psvd=30％

圖3和圖4分別顯示了SNR為150和10時，對不同駐留函數(shù)和不同血流量值，選擇不同的閾值奇異值分解法對腦血流量的估計結(jié)果。對不同駐留函數(shù)不同血流量值所選閾值如表1所示。

對于圖3的高信噪比情況，與圖1的結(jié)果相比，可發(fā)現(xiàn)，在高血流量時，奇異值分解依然能夠很好估計腦血流量。這說明在信噪比高時，奇異值分解法可以有效地估計腦血流量，血流量大時的標準差明顯大于血流量小時。

對于圖4的低信噪比情況，與圖2的結(jié)果相比，可發(fā)現(xiàn)不管是在高血流量和低血流量時，估計結(jié)果均得到明顯改善。說明信噪比小時，奇異值分解法對腦血流量的估計依然有效，但標準差隨著血流量的增大而增大。

可見，無論信噪比高低，奇異值分解法都可以很好地估計腦血流量。在信噪比較高時，估計結(jié)果的標準差小于低信噪比時的標準差；信噪比相同時，低血流量時估計結(jié)果的標準差小于高血流量時的標準差。這與?stergaard L等的仿真結(jié)果相符［4］。

表1 對不同駐留函數(shù)不同血流量值所選閾值（％）Tab.1 The threshold for different residue function and different CBF

圖3 SNR=150時選擇不同閾值的腦血流量估計結(jié)果。（a）腦血流量估計值；（b）標準差Fig.3 Estimated CBF values with SNR=150 and differentthreshold.（a）estimated CBF values；（b）standard deviation

3.3 動脈輸入函數(shù)（AIF）延遲和失真對腦血流量估計的影響

圖4 SNR=10時選擇不同閾值的腦血流量估計結(jié)果。（a）腦血流量估計值；（b）標準差Fig.4 Estimated CBF values with SNR=10 anddifferentthreshold.（a）estimated CBF values；（b）standard deviation

圖5為AIF延遲時奇異值分解法對腦血流量的估計結(jié)果。結(jié)果顯示，對于指數(shù)型和三角型駐留函數(shù)，當(dāng)血流量低于30 mL/100 mg/min時，AIF延遲不會對腦血流量估計結(jié)果造成影響，當(dāng)血流量高于30mL/100 mg/min時，AIF延遲會造成腦血流量的明顯低估，且血流量越大，低估情形越嚴重；對于盒型駐留函數(shù)，低估情形不明顯。這是因為盒型駐留函數(shù)過于理想，不能很好地描述示蹤劑到達組織后隨著時間還留在組織內(nèi)的幾率。

圖6為AIF失真時奇異值分解法對腦血流量的估計結(jié)果。結(jié)果顯示AIF失真對估計結(jié)果的影響并不明顯。

圖5 AIF延遲時的估計結(jié)果Fig.5 Estimated CBF values with AIF delay

圖6 AIF失真時的估計結(jié)果Fig.6 Estimated CBF values with AIF distortion

3.4 延遲效應(yīng)的修正

圖7為修正延遲效應(yīng)前后腦血流量估計結(jié)果的比較，所用駐留函數(shù)為指數(shù)型?？梢钥闯觯貉舆t時間的長短對腦血流量低估程度影響不大；腦血流量越大時，低估程度越嚴重，當(dāng)血流量為70 mL/100mg/min時，腦血流量會被低估30％左右，當(dāng)血流量為40 mL/100 mg/min時，腦血流量會被低估20％左右；修正延遲效應(yīng)后低估程度減小到±5％，得以改善，但是方差比修正延遲效應(yīng)前明顯增大，特別在延遲時間較短時。若駐留函數(shù)為三角型，修正延遲效應(yīng)前后仿真實驗結(jié)果與指數(shù)型駐留函數(shù)相似。

圖7 修正延遲效應(yīng)前后估計結(jié)果的比較。（a）CBF=70 mL/100 mg/min；（b）CBF=60 mL/100 mg/min；（c）CBF=50 mL/100 mg/min；（d）CBF=40 mL/100 mg/minFig.7 Results of simulations for comparison between delay-corrected CBF and uncorrected CBF.（a）CBF=70 mL/100 mg/min；（b）CBF=60 mL/100 mg/min；（c）CBF=50 mL/100 mg/min；（d）CBF=40 mL/100 mg/min

4 結(jié)論

本研究對奇異值分解法估計腦血流量進行了理論推導(dǎo)，采用了指數(shù)型、三角型和盒型三種駐留函數(shù)進行了仿真實驗，分別對不同的信噪比、不同的閾值選擇、動脈輸入函數(shù)的延遲與失真進行了研究。

仿真實驗表明，對3種不同的駐留函數(shù)，不同的血流量采用相同的閾值會造成部分血流量估計結(jié)果不準確。這是因為血流量值大時的示蹤劑濃度-時間曲線的峰值較大，選擇相同的閾值會造成高血流量時有用信息的濾除或低血流量時噪聲不能完全濾除。對不同的血流量值采用不同的閾值可達到更好的效果。對3種駐留函數(shù)，不同的血流量采用不同的閾值，無論信噪比高低，奇異值分解法都可以有效地估計腦血流量。仿真研究表明，動脈輸入函數(shù)發(fā)生失真時，該方法依然可以有效地估計腦血流量，即該方法對AIF失真不敏感。但是，該方法對動脈輸入函數(shù)發(fā)生延遲比較敏感，尤其是當(dāng)腦血流量較高時，該方法會造成腦血流量的低估。研究發(fā)現(xiàn)，低估的程度與血流量的大小有關(guān)，血流量越大低估情形越嚴重，而低估的程度與延遲時間的長短沒有關(guān)系。經(jīng)過對延遲效應(yīng)影響的修正，低估情形得以明顯改善，但方差有所增加，特別是在延遲時間較短時。

本研究的仿真結(jié)果表明，通過不同閾值的選擇和對延遲時間的修正，奇異值分解法可對MR灌注成像腦血流量進行有效地估計，為下一步的臨床數(shù)據(jù)處理打下基礎(chǔ)，為腦部疾病的臨床診斷及治療效果的評估提供手段。

［1］Petrella JR，Provenzale JM.MR perfusion imaging of the brain：

techniques and applications ［J］.American Roentgen Ray Society，2000，175：207-219.

［2］Wong JC，Provenzale JM，Pertrella JR.Perfusion MR imaging of brain neoplasms［J］.American Roentgen Ray Society，2000，174：1147-1157

［3］俎棟林.核磁共振成像學(xué)［M］.北京：高等教育出版社，2004

［4］?stergaard L， WeisskoffRM， CheslerDA， etal.High resolution measurement of cerebral blood flow using intravascular tracer boluspassages.PartI：Mathematicalapproach and statistical analysis［J］.Magn Reson Med，1996，36（5）：715-725.

［5］Wirestam R，Andersson L，?stergaard L，et al.Assessment of regional cerebral blood flow by dynamic susceptibility contrast MRI using different deconvolution techniques［J］.Magn Reson Med，2000，43（5）：691-700.

［6］Calamante F，Gadian DG，Connelly A.Delay and dispersion effects in dynamic susceptibility contrast MRI：simulations using singular value decomposition ［J］.Magn Reson Med，2000，44（3）：466-473.

［7］Wu O，?stergaard L，Koroshetz WJ，et al.Effects of tracer arrival time on flow estimates in MR perfusion-weighted imaging［J］.Magn Reson Med，2003，50（4）：856-864.

［8］Smith MR，Lu H，F(xiàn)rayne R.Signal-to-noise ratio effects in quantitative cerebralperfusion using dynamic susceptibility contrast agents［J］.Magn Reson Med，2003，49（1）：122-128.

［9］Smith MR，Lu H，Trochet S，et al.Removing the effect of SVD algorithmic artifacts present in quantitative MR perfusion studies［J］.Magn Reson Med，2004，51（4）：631-634.

［10］李承鋒.利用曲線擬合找到血流灌注磁振影像的血液動力參數(shù)［D］.臺灣：國立陽明大學(xué)，2006.

［11］Weisskoff RM，Zuo CS，Boxerman JL，et al.Micro-scopic susceptibility variation and transverse relaxation：theory and experiment［J］.Magn Reson Med，1994，31：601-610.

［12］Fisel CR，Ackerman JL，Buxton RB，et al.MR contrast due to microscopically heterogeneous magnetic susceptibility：numerical simulations and applications to cerebral physiology ［J］.Magn Reson Med，1991，17：348-356.

［13］Unal S，Rohit S.Cerebral blood flow estimation from perfusionweighted MRI using FT-based MMSE filtering method ［J］.Magnetic Resonance Imaging，2008，26：313-322.