基于小波包變換的特征提取方法分析fMRI數(shù)據(jù)

支聯(lián)合 譚素敏 楊建國

1（周口師范學(xué)院物理與電子工程系，周口 466001）

2（周口師范學(xué)院附屬醫(yī)院，周口 466001）

3（哈爾濱工業(yè)大學(xué)威海校區(qū)汽車工程學(xué)院，威海 264209）

引言

由于具有無創(chuàng)、可重復(fù)和時(shí)空分辨率高等特點(diǎn)，功能磁共振成像技術(shù)（functional magnetic resonance imaging，fMRI）在人腦解剖定位、功能成像和網(wǎng)絡(luò)連接研究中起著日益關(guān)鍵的作用［1－2］。但是，fMRI數(shù)據(jù)里除包含激活信號外，還混雜有高頻噪聲、低頻漂移等干擾。因此，設(shè)計(jì)靈敏度高、特異度好的fMRI數(shù)據(jù)處理方法成為重要的研究課題。其中，靈敏度是識別激活體素的能力，取決于去除干擾能力的強(qiáng)弱；特異度是區(qū)分激活體素和非激活體素的能力，決定于去除干擾時(shí)對激活信號的保護(hù)程度。

統(tǒng)計(jì)參數(shù)圖（statistical parametric mapping，SPM）是目前國際權(quán)威的fMRI分析方法［3］，用試驗(yàn)刺激函數(shù)與血流動力學(xué)函數(shù)卷積模擬試驗(yàn)刺激誘導(dǎo)的激活信號，用高通濾波去除低頻漂移，再基于廣義線性模型對體素?cái)?shù)據(jù)進(jìn)行分析，最新版本為SPM8。具有多尺度分析功能的小波變換也廣泛應(yīng)用于fMRI數(shù)據(jù)在內(nèi)的醫(yī)學(xué)信號處理［4－6］。比如，基于小波去噪理論［7］去除fMRI數(shù)據(jù)的高頻噪聲［8］，基于多尺度分析特性去除fMRI數(shù)據(jù)的低頻漂移［9］，在小波域?qū)MRI數(shù)據(jù)進(jìn)行廣義線性模型建模［10］等。

與上述方法僅去除單一高頻或低頻干擾不同，多尺度特征提取（multiscale feature extraction，MFE）基于離散小波變換（discrete wavelet transform，DWT）分割fMRI數(shù)據(jù)的頻譜，進(jìn)而通過重構(gòu)激活信號存在的子帶同時(shí)去除高頻和低頻干擾［11］，顯著提高了后續(xù)統(tǒng)計(jì)檢測性能。理論上，小波包變換（wavelet packet transform，WPT）比小波變換對頻帶的分割更精細(xì)，因此更適合MFE提取激活信號和去除干擾。為此，基于WPT構(gòu)建新MFE，并設(shè)計(jì)矩陣算法代替逐體素迭代算法快速提取激活信號，同時(shí)與原有MFE和SPM8進(jìn)行比較。

1 材料和方法

1.1 聽覺fMRI試驗(yàn)

一個(gè)被試的聽覺fMRI試驗(yàn)數(shù)據(jù)下載于SPM網(wǎng)站（www.fil.ion.ucl.ac.uk/spm），并被許可使用。試驗(yàn)為組塊型設(shè)計(jì)，任務(wù)組塊刺激是給被試聽雙音節(jié)單詞，語速為60個(gè)/min。試驗(yàn)以靜息組塊開始，任務(wù)與靜息兩種組塊交替進(jìn)行，每個(gè)組塊持續(xù)42 s。用2 Tesla Siemens掃描機(jī)以f0=1/7Hz的抽樣頻率獲取一個(gè)被試共96幀全腦BOLD/EPI圖像。每幀圖像大小為64體素×64體素×64體素，體素大小為3 mm×3 mm×3 mm。前12幀圖像因磁場穩(wěn)定性原因舍棄不用，剩余84幀圖像則進(jìn)行SPM8方法常規(guī)預(yù)處理，包括被試運(yùn)動矯正、標(biāo)準(zhǔn)化到 M NI解剖空間、空間高斯平滑及去除腦外體素。

試驗(yàn)刺激的周期為 T =84 s，時(shí)寬為42 s，因此激活信號的頻率組分除直流外，還應(yīng)包括基頻1/T=0.011 9Hz和三倍頻 3 /T=0.035 7Hz、五倍頻5/T=0.059 5Hz等。圖1（a）顯示了SPM8模擬的激活信號X的頻譜。

1.2 方法

1.2.1 小波包變換

WPT基于DWT發(fā)展而來，包括分解和重構(gòu)兩個(gè)互逆的運(yùn)算，為易于理解先敘述DWT的原理。對于確定的小波低通濾波器L、高通濾波器H和下采樣算子D0，一個(gè)體素時(shí)間序列數(shù)據(jù)y經(jīng)J級DWT分解后變換為分布在J+1個(gè)時(shí)頻尺度里的小波系數(shù)集合cy=［aJ，dJ，dJ－1，…，d1］

式中，分解深度j=1，2，…，J；aj和dj分別為第j級低頻尺度和高頻尺度。與DWT分解只分割低頻尺度不同，WPT對高頻尺度和低頻尺度同時(shí)進(jìn)行分割

式中，分解深度j=1，2，…，J；在任一深度j，i=1，2，…，I－1，I， 其中 I =2j－1；sij－1表示第j－1層的第i個(gè)小波包尺度。上述分解都稱為多尺度分析，把體素頻帶分割成了一系列時(shí)頻特性各異的尺度（子帶）。對于理想小波濾波器，aJ的頻帶為［0，f0/2J+1），dj的為（f0/2j+1，f0/2j］，WPT的 2J個(gè)尺度則平均分割體素頻帶［0，f0/2］?？梢?，分解同樣的深度時(shí)WPT對頻帶的分割更為精細(xì)，這是選用WPT代替DWT構(gòu)建新MFE的原因。但因?yàn)榉指罹?xì)，WPT耗時(shí)因此更多。

重構(gòu)是分解的逆過程，DWT的重構(gòu)為

式中，U為隔點(diǎn)插零算子，和分別為小波低通和高通重構(gòu)濾波器，這時(shí)的j=J，J－1，…，1。而WPT則為

上述是對全部尺度的重構(gòu)，如果重構(gòu)時(shí)保留特定尺度里的系數(shù)不變，而將剩余尺度里的系數(shù)全部置零則稱為特征尺度重構(gòu)。這種重構(gòu)雖則多了系數(shù)置零的步驟，但能夠有選擇地從y中提取感興趣的子帶成分。

1.2.2 多尺度特征提取fMRI數(shù)據(jù)

WPT更精細(xì)的頻帶分割能力為利用特征尺度重構(gòu)提取fMRI數(shù)據(jù)里的激活信號提供了更好的技術(shù)支持，新MFE包括3個(gè)步驟。

步驟1：識別特征尺度。激活信號存在的尺度為特征尺度，其它尺度為干擾尺度。取J=4用WPT把模擬的激活信號X分解在16個(gè)尺度里，并逐尺度地進(jìn)行重構(gòu)。對比每個(gè)重構(gòu)尺度的頻譜與X的頻譜，發(fā)現(xiàn)X的成分主要分布在除、和外的尺度里。這樣，就把上述3個(gè)尺度及確定為干擾尺度，剩余的12個(gè)為特征尺度。盡管尺度包含X的直流成分，但體素?cái)?shù)據(jù)的這一尺度里也包含低頻漂移，故也視為干擾尺度。對X進(jìn)行特征尺度重構(gòu)，得到的信號為E1X，其頻譜顯示在圖1（b）里。

步驟2：提取 y 里包含的激活信號。也即對 y進(jìn)行特征尺度重構(gòu)，得到的信號為E1y。

圖1 3種方法分析聽覺數(shù)據(jù)涉及到的頻譜。（a）SPM8；（b）E1X；（c）尺度d31；（d）E2XFig.1 The frequency spectrums used by three methods for the auditory data.（a）SPM8；（b）E1X；（c）the scale d31；（d）E2X

fMRI數(shù)據(jù)里體素具有海量性質(zhì)，WPT速度又慢，因此常規(guī)基于迭代算法的一個(gè)體素接著一個(gè)體素地進(jìn)行特征尺度重構(gòu)非常耗時(shí)。為此，設(shè)計(jì)基于WPT的矩陣算法一次性地對諸多體素進(jìn)行特征尺度重構(gòu)式中，Y是由體素?cái)?shù)據(jù)作為每一列組成的數(shù)據(jù)矩陣，特征提取矩陣M則通過對N階單位方陣的每一列都進(jìn)行特征尺度重構(gòu)得到，N為體素?cái)?shù)據(jù)長度，在本研究為84；eY為提取的激活信號矩陣，每一列對應(yīng)一個(gè)體素。

步驟3：檢測激活體素。對E1X和eY里的每一列數(shù)據(jù)進(jìn)行相關(guān)分析，并將相關(guān)系數(shù)r用Fisher'sZ變換轉(zhuǎn)化為Z值

對于基于DWT的MFE方法，第一步也是識別特征尺度。同樣取J=4，分析頻譜發(fā)現(xiàn)X的成分主要分布在尺度d31（即d3和d1的組合）里，但重構(gòu)這兩個(gè)尺度得到的頻譜顯示了嚴(yán)重的混頻，即出現(xiàn)了X沒有的頻率成分，見圖1（c）。這是把尺度 a4和d42里的系數(shù)設(shè)置為零的緣故。為避免混頻，取J=5且只把a(bǔ)5視為干擾尺度，而把d54321視為特征尺度，得到的重構(gòu)信號記為E2X，其頻譜顯示在圖1（d）里。確定特征尺度之后，剩余的兩個(gè)步驟與WPT的相同。兩種方法的程序編寫都基于Matlab7.2小波工具箱，小波都用sym2。

對于SPM8，采用其網(wǎng)站提供的方法進(jìn)行分析，用到的高通濾波的截止頻率為其缺省的1/128Hz，基于廣義線性模型用t-檢驗(yàn)逐體素地檢驗(yàn)任務(wù)與靜息兩狀態(tài)的差異，最后的t-值通過該軟件自帶的程序轉(zhuǎn)化為Z值。

選取一定的顯著性水平作為閾值，對3種方法得到的Z值進(jìn)行閾值化。因多重統(tǒng)計(jì)檢驗(yàn)的原因，作為閾值的顯著性水平都基于高斯隨機(jī)場理論進(jìn)行修正，修正由SPM8完成。與聽覺有關(guān)的軸向坐標(biāo)z=6 mm、10 mm和14 mm腦層上的激活體素按放射學(xué)約定的方式顯示在MNI解剖圖上。

2 結(jié)果

圖2顯示了基于WPT的新MFE、原有的基于DWT的MFE和SPM8分析聽覺fMRI試驗(yàn)數(shù)據(jù)3層腦區(qū)的激活Z-圖，圖中最上面3行選取的閾值為顯著性水平P=0.001，最下面一行P=0.05，它們都經(jīng)過修正?？梢钥闯?，盡管3種方法的分析機(jī)理不同，但檢測到的激活區(qū)的位置都相同，并且都在公認(rèn)的聽覺區(qū)，因此3種方法的特異度相同，分析結(jié)果是可靠的。但同時(shí)新MFE檢測到的激活區(qū)最大，原有MFE的次之，SPM8的最小。因此，新MFE的分析結(jié)果在3種方法中是最好的。

表1更詳細(xì)地提供了這3種方法檢測圖2中3層腦數(shù)據(jù)區(qū)得到的激活體素個(gè)數(shù)。與圖2一樣，表1基于WPT和DWT的MFE的閾值都選為修正的P=0.001，SPM8的為修正的P=0.05。顯然，對于每層數(shù)據(jù)而言，基于WPT的新MFE檢測出的激活體素個(gè)數(shù)都是最多的。進(jìn)一步計(jì)算可知，基于WPT的新MFE檢測的激活體素個(gè)數(shù)比原有MFE多13.2％，比SPM8多30.8％。

因此，綜合上述分析結(jié)果可知，基于WPT的MFE的分析性能在3種方法中是最優(yōu)的。

表1 3種方法分析聽覺數(shù)據(jù)檢測的激活體素個(gè)數(shù)Tab.1 Numbers of active pixels detected by three methods for auditory data

圖2 3種方法分析聽覺fMRI試驗(yàn)一個(gè)被試3層數(shù)據(jù)的激活Z-圖（圖中顯著性水平分別為P=0.001（上三行）和P=0.05（最下一行），這些值都經(jīng)過高斯隨機(jī)場理論修正）Fig.2 Theactive Z-mapsgenerated bythethree methods for the auditory fMRI data of a subject（The significant level are P=0.001（the first three rows）and P=0.05（the bottom row）respectively，which are corrected based on Gaussian random field）

3 討論

WPT顯著提高M(jìn)FE分析性能的根本原因，在于WPT的時(shí)頻尺度在多尺度分析方法中是最窄的。MFE分析體素信號的根本是借助時(shí)頻特性各異的尺度把信號頻帶分割成一系列的子帶，在重構(gòu)時(shí)舍棄干擾尺度對應(yīng)的子帶和保留特征尺度對應(yīng)的子帶。顯然，如要獲得較高的靈敏度，特征尺度里必須包含較少的干擾成分，否則重構(gòu)時(shí)這些干擾無法去除；同樣，為了獲得較高的特異度，干擾尺度也應(yīng)該包含較少的激活信號成分，因?yàn)檫@些成分將在重構(gòu)時(shí)舍棄。但在fMRI數(shù)據(jù)的頻譜里，激活信號的頻率組分不是全部都集中在一個(gè)區(qū)段里，而是離散地分布于各干擾成分之間。這樣，為了達(dá)到上述目的，MFE采用的時(shí)頻尺度的頻寬必須盡可能地窄，以便在精確捕捉激活信號頻率組分的同時(shí)包含較少的干擾成分。目前常用的多尺度分析方法，除了本文涉及到的WPT和 D WT，還包括提升小波變換和平穩(wěn)小波變換，但尺度最窄的是 W PT。因此，WPT賦予MFE更好的分析性能。

例如，原有MFE方法基于DWT取J=5把體素頻帶分割為6個(gè)尺度，其中d1的頻帶范圍為［f0/4，f0/2］，涵蓋體素頻譜的一半，里面雖有激活信號的五倍頻成分，但同時(shí)包含大量的高頻噪聲。這樣，如把d1選為干擾尺度，則會丟失激活信號，從而損害特異度；如選為特征尺度，則無法去除高頻噪聲，因此損害靈敏度。所以，對于原有 M FE來說，無論是把d1確定為特征尺度還是干擾尺度都不是最優(yōu)選擇。相反，新MFE基于WPT取J=4時(shí)，就已把上述頻帶范圍又細(xì)分割成了從一直到共8個(gè)尺度，而且尺度、和里沒有激活信號成分。于是，通過把這3個(gè)尺度選為干擾尺度新MFE就獲得了比原有MFE和SPM8更高的靈敏度。而且，由于這樣做沒有損害激活信號的五倍頻成分，所以新MFE又同時(shí)維持了與另外兩種方法一樣的特異度。

除了分析性能方面的考慮外，消耗時(shí)間少的方法顯然更容易為使用者接受。由于WPT比DWT分割出的尺度多，因此運(yùn)算時(shí)間相對要多，這對于分析海量性質(zhì)的fMRI數(shù)據(jù)來說是個(gè)缺點(diǎn)。對于本文的聽覺數(shù)據(jù)，在同樣配置的計(jì)算機(jī)條件下，基于WPT和DWT矩陣算法的MFE所用時(shí)間都為31 s。SPM8分析fMRI數(shù)據(jù)的機(jī)理與MFE不同，但消耗在超參數(shù)估計(jì)、參數(shù)估計(jì)和平滑估計(jì)上的時(shí)間也為77 s。相反，采用逐體素的方法反復(fù)迭代WPT分解、干擾尺度系數(shù)置零以及WPT重構(gòu)這3個(gè)步驟，總共消耗了2 918 s，約合48.5 min，在實(shí)際中這顯然是不能接受的。因此，本文給出的矩陣算法保證了新MFE方法在實(shí)際應(yīng)用中分析海量數(shù)據(jù)的可行性。

和DWT一樣，WPT對小波的選取比較敏感，不同的小波可能給出不同的結(jié)果。比如，混頻是基于DWT和WPT的特征提取方法的一個(gè)共性問題，源于小波濾波器的非理想特性、下抽樣不滿足抽樣定理等，就與小波的選取有關(guān)［12］。文獻(xiàn)［12］研究發(fā)現(xiàn)，小波濾波器的系數(shù)越多，混頻現(xiàn)象就越不明顯。所以，對于一個(gè)具體的實(shí)際數(shù)據(jù)而言，合適小波的選取涉及到分析性能、時(shí)間消耗和使用者對小波知識的了解等諸多因素，需要綜合考量，因此是一個(gè)復(fù)雜的問題。在多尺度分析中，只有平穩(wěn)小波變換對小波的選取不太敏感，采用這一變換可以免去小波選取的環(huán)節(jié)，但這一變換方法與DWT一樣對頻譜的分割沒有WPT精細(xì)。所以，未來的研究方向是設(shè)計(jì)一種既是平穩(wěn)的又是小波包的多尺度分析方法，即平穩(wěn)小波包變換，據(jù)此獲得一種既能精細(xì)分割頻譜又對小波的選取不太靈敏的fMRI數(shù)據(jù)分析方法。

4 結(jié)論

總之，在目前的多尺度分析技術(shù)范圍內(nèi)，WPT賦予MFE更好的靈敏度和特異度來分析fMRI數(shù)據(jù)，本研究引入的WPT矩陣算法較逐體素的反復(fù)迭代算法顯著提高了MFE的運(yùn)算速度，所提出的是一種更為實(shí)用的fMRI數(shù)據(jù)分析方法。

［1］堯德中，羅程，雷旭，等.腦成像與腦連接［J］.中國生物醫(yī)學(xué)工程學(xué)報(bào)，2011，30（1）：6－10.

［2］包尚聯(lián).腦功能成像物理學(xué)［M］.鄭州：鄭州大學(xué)出版社，2006：200－207.

［3］Friston KJ，HolmesAP，Worsley KJ，etal. Statistical parametric maps in functional imaging：a general linear approach［J］.Hum Brain Mapping，1995，2（4）：189 －210.

［4］Bullmore E，F(xiàn)adili J，Maxim V，et al.Wavelets and functional magnetic resonance imaging ofthe human brain ［J］.NeuroImage，2004，23（1）：S234－S249.

［5］ 代 煜，張建勛，雪原.基于小波變換的脊柱振動特征分析［J］.中國生物醫(yī)學(xué)工程學(xué)報(bào)，2012，31（3）：461－465.

［6］ 王 永軒，邱天爽，劉蓉.基于小波分析方法的腦電誘發(fā)電位單導(dǎo)少次提?。跩］.中國生物醫(yī)學(xué)工程學(xué)報(bào)，2011，30（1）：34－39.

［7］ D onoho DL，Johnstone IM.Adapting to unknown smoothness via wavelet shrinkage ［J］.J Am Stat Assoc，1995，90（432）：1200－1224.

［8］ W ink AM，Roerdink JBTM.Denoising functional MR images：a comparison of wavelet denoising and Gaussian smoothing ［J］.IEEE Trans Med Imaging，2004，23（3）：374－387.

［9］ M eyerFG. W avelet-based estimation ofa semiparametric generalized linear model of fMRI time-series［J］.IEEE Trans Med Imaging，2003，22（3）：315－322.

［10］ V an De Ville D，Seghier ML，Lazeyras F，et al.WSPM：wavelet-based statistical parametric mapping ［ J］.NeuroImage，2007，37（4）：1205－1217.

［11］ 支 聯(lián)合，李玉曉，趙書俊，等.基于離散小波變換的 f MRI數(shù)據(jù)特征提?。跩］.中國醫(yī)學(xué)影像技術(shù)，2010，26（6）：1151－1154.

［12］ 楊 建國.小波分析及其工程應(yīng)用［M］.北京：機(jī)械工業(yè)出版社，2005：107－128.