近似熵與樣本熵在靜息態(tài)fMRI數(shù)據(jù)中的對比研究

陳 慧， 孟攀婷， 朱登輝， 趙 晨， 張亮亮

(1.鄭州輕工業(yè)大學(xué) 工程訓(xùn)練中心， 河南 鄭州 450000; 2.鄭州輕工業(yè)大學(xué) 計算機(jī)與通信工程學(xué)院， 河南 鄭州 450000)

0 引言

隨著功能磁共振成像(functional magnetic resonance imaging,fMRI)的迅速發(fā)展，越來越多的人開始進(jìn)行fMRI的研究.大腦是一個非常復(fù)雜的非線性系統(tǒng)，近幾年來越來越多的非線性方法被廣泛應(yīng)用于fMRI中以探索大腦功能的非線性特征信息，其中用之較多的是通過熵[1-4]從時間序列的復(fù)雜性角度來探索大腦的功能活動.

Sokunbi等[5-7]的研究表明將近似熵(approximate entropy,ApEn)方法應(yīng)用于fMRI中是可行的、有效的.Sokunbi等隨后的一系列研究嘗試將樣本熵(sample entropy,SampEn)方法應(yīng)用于fMRI中，結(jié)果發(fā)現(xiàn)SampEn方法在區(qū)分正常人與疾病(注意力不集中癥、精神分裂癥)患者、年輕人與老年人方面均具有較好的效果.以上研究表明ApEn方法與SampEn方法在fMRI研究中均具有較好的應(yīng)用前景，然而究竟采用哪種算法分析fMRI數(shù)據(jù)更有效，需要進(jìn)一步的討論與研究.

基于此，本文針對腸易激綜合癥(irritable bowel syndrome,IBS)與健康對照(healthy control,HC)兩組人的靜息態(tài)fMRI數(shù)據(jù)，采用ApEn方法與SampEn方法進(jìn)行分析，計算不同參數(shù)r下的ApEn與SampEn的均值差距率并進(jìn)行對比，探討哪種算法更適用于fMRI數(shù)據(jù)的分析.

1 近似熵與樣本熵方法

1.1 影像學(xué)fMRI數(shù)據(jù)參數(shù)

影像學(xué)數(shù)據(jù)是在被試處于靜息閉眼狀態(tài)(未受到任何刺激，也不執(zhí)行任何任務(wù))下，通過德國西門子公司的3T磁共振掃描儀采用梯度平面回波成像序列得到的，其參數(shù)如下，層厚度5 mm、矩陣大小64×64、回波時間30 ms、重復(fù)時間2 s、翻轉(zhuǎn)角度90 °、層內(nèi)分辨率3.75×3.75 mm2.每個被試的掃描過程需要6 min，共采集180個時間點.在掃描過程中，被試被要求保持清醒并放松的狀態(tài)，同時使用烏龍頭線圈和泡沫墊子固定被試頭部避免頭動.在本研究中，影像學(xué)數(shù)據(jù)包含16例IBS患者與年齡、性別相匹配的16例HC，其統(tǒng)計學(xué)信息如表1所示.

表1 兩組被試的人口統(tǒng)計學(xué)信息

1.2 影像學(xué)的數(shù)據(jù)預(yù)處理

本文對影像學(xué)數(shù)據(jù)的預(yù)處理是基于中科院心理研究所開發(fā)的DPABI[8]進(jìn)行的，主要步驟包括：①將fMRI的原始DICOM數(shù)據(jù)轉(zhuǎn)換為NIFTI數(shù)據(jù)；②去除體素的時間序列的前10個時間點，保留后170個時間點，避免剛開始掃描儀的不穩(wěn)定及被試的不適應(yīng)性等因素對結(jié)果的影響；③進(jìn)行時間層校正，消除層與層的掃描時間的差異；④進(jìn)行頭動校正，去除頭動移動大于1.5 mm或者旋轉(zhuǎn)角度大于1.5 °的被試數(shù)據(jù)，避免在掃描期間頭動對結(jié)果的影響；⑤進(jìn)行空間標(biāo)準(zhǔn)化,消除不同個體大腦形態(tài)的差異；⑥進(jìn)行去線性漂移及帶通濾波，去除掃描儀及生理活動產(chǎn)生的噪聲，提高信噪比，其中帶通濾波頻率為0.01～0.08 Hz.

1.3 近似熵方法在fMRI中的實現(xiàn)過程

本研究中近似熵的計算是基于MATLAB平臺和SPM工具包實現(xiàn)的，經(jīng)過預(yù)處理后的個體腦圖像均包含70 831個3.75×3.75 mm3大小的體素，每一個體素對應(yīng)著一個時間序列，對體素逐個進(jìn)行ApEn求解，具體步驟如下：

①設(shè)預(yù)處理后fMRI個體數(shù)據(jù)的腦圖像有n個(即時間序列的長度為n)，全腦有l(wèi)個體素；

②提取全腦體素的時間序列，矩陣形式如下：

其中，每一列表示一個體素的時間序列.

③針對第t個體素的時間序列X(t)=X(:,t)，將其重構(gòu)成m=2維向量，向量形式如下：

Xi(t)=(xi(t),xi+1(t),…,xi+m-1(t))T,

i=1,2,…,n-m+1

④計算向量Xi(t),i=1,2,…,n-m+1與所有向量Xj(t),j=1,2,…,n-m+1之間的距離，此時的距離為兩向量對應(yīng)元素差值絕對值的最大值，其計算公式如下：

dij=max|xi+k(t)-xj+k(t)|,

k=0,1,…,m-1

(1)

(2)

⑦維數(shù)m增加1，重復(fù)步驟③～⑥，求得Φm+1(t)；

⑧將求得的Φm(t)和Φm+1(t)作差，其結(jié)果即為ApEn(t)，其計算公式如下：

ApEn(t)=Φm(t)-Φm+1(t)

(3)

1.4 樣本熵方法在fMRI中的實現(xiàn)過程

本研究中樣本熵的計算是基于MATLAB平臺和SPM工具包實現(xiàn)的.針對預(yù)處理后的個體腦圖像，對體素逐個進(jìn)行SampEn求解，具體步驟如下：

①設(shè)預(yù)處理后fMRI個體數(shù)據(jù)的腦圖像有n個(即時間序列的長度為n)，全腦有l(wèi)個體素；

②提取全腦體素的時間序列，矩陣形式如下：

其中，每一列表示一個體素的時間序列.

Xi(t)=(xi(t),xi+1(t),…,xi+m-1(t))T,

i=1,2,…,n-m+1

④計算向量Xi(t),i=1,2,…,n-m與其他向量Xj(t),j=1,2,…,n-m之間的距離，此時的距離為兩向量對應(yīng)元素差值絕對值的最大值，其計算公式如下：

dij=max|xi+k(t)-xj+k(t)|,

k=0,1,…,m-1,i≠j

(4)

⑥計算所有向量中小于閾值的數(shù)目的平均值，記為Φm(t)，其計算公式如下：

水土流失、植被破壞是我國泥石流多發(fā)的原因之一。在泥石流等地質(zhì)災(zāi)害高發(fā)區(qū)，應(yīng)加大環(huán)境保護(hù)力度，實行治水、治土和工程防護(hù)的綜合治理策略，是為治本之策。我國環(huán)境綜合治理取得了較為顯著的效果，應(yīng)進(jìn)一步推進(jìn)。目前主要的幾種綜合治理模式如表1所示。

(5)

⑦維數(shù)m增加1，重復(fù)步驟③～⑥，求得Φm+1(t)；

⑧將求得的Φm(t)和Φm+1(t)取對數(shù)并作差，其結(jié)果即為SampEn(t)，其計算公式如下：

SampEn(t)=lnΦm(t)-lnΦm+1(t)

(6)

1.5 均值差距率

熵用于衡量時間序列的信息發(fā)生率，其值的大小反映了系統(tǒng)的復(fù)雜性程度.已有研究表明，年齡和疾病特征體現(xiàn)在復(fù)雜生理系統(tǒng)的動態(tài)復(fù)雜度缺失，即熵的降低.因此，IBS患者組與HC組兩組人的熵必然會存在差異.本文通過均值差距率v來量化IBS患者組與HC組兩組人熵的差異，其值的大小可反映了HC組與IBS組兩組人熵的差異程度，其值越大，HC組與IBS組兩組人熵的差異越明顯，其計算公式如下：

(7)

式(7)中：mHC表示HC組ApEn或SampEn的平均值，mIBS表示IBS組ApEn或SampEn的平均值.

1.6 fMRI數(shù)據(jù)處理及統(tǒng)計分析

采用步驟1.3的ApEn方法與步驟1.4的SampEn方法計算不同參數(shù)r下的兩組被試所有體素的ApEn與SampEn，并計算所有被試ApEn與SampEn的全腦平均，每個被試對應(yīng)一個ApEn與SampEn的全腦平均.采用雙樣本t檢驗對HC組與IBS組兩組被試的全腦平均ApEn與SampEn進(jìn)行統(tǒng)計分析(本文認(rèn)為當(dāng)p<0.05，兩組被試的全腦平均ApEn值與全腦SampEn值存在顯著性差異).采用式(7)計算HC組與IBS組的ApEn與SampEn的均值差距率.

為了使HC組與IBS組的均值差距率更加明顯，本文同時進(jìn)行了如下步驟：①對不同參數(shù)r下的兩組被試所有體素的ApEn與SampEn進(jìn)行雙樣本t檢驗(p<0.05)，篩選出具有顯著性差異的體素；②計算具有顯著性差異的體素的平均ApEn與平均SampEn；③采用式(7)計算HC組與IBS組的ApEn與SampEn的均值差距率.

根據(jù)上述計算結(jié)果，可以得到ApEn方法與SampEn方法最優(yōu)的參數(shù)r.計算將該最優(yōu)參數(shù)r下的所有體素的ApEn與SampEn，生成相應(yīng)的ApEn圖與SampEn圖并進(jìn)行平滑處理.采用雙樣本t檢驗(GRF校正(p<0.05,差異簇大小>20))進(jìn)行統(tǒng)計分析，確定兩組被試的差異性區(qū)域，以驗證ApEn方法與SampEn方法的合理性.

2 結(jié)果與討論

2.1 近似熵與樣本熵的均值差距率分析

本文采用ApEn方法與SampEn方法計算不同參數(shù)r下的所有被試ApEn與SampEn的全腦平均，其值隨r的變化曲線如圖1所示.采用雙樣本t檢驗對不同參數(shù)r下的全腦平均ApEn與SampEn進(jìn)行統(tǒng)計分析，其結(jié)果如表2所示.從表2可以看出，ApEn是在r=0.2～0.4下呈現(xiàn)顯著性組差異，SampEn是在r=0.15～0.55下呈現(xiàn)出顯著性組差異(需要指出，在r=0.1，由于選取r值過小，使得該參數(shù)r下存在體素的SampEn值為Inf的情況，顯然這是不符合實際的).

圖1 ApEn與SampEn的全腦平均隨r的變化曲線

表2 兩組被試的全腦平均ApEn與SampEn在不同參數(shù)r下的統(tǒng)計學(xué)差異

采用均值差距率計算公式分別計算不同參數(shù)r下的ApEn與SampEn的均值差距率(未經(jīng)過t檢驗篩選差異體素)，其值隨r的變化曲線如圖2所示.為了使HC組與IBS組的均值差距率更加明顯，更好地評估ApEn和SampEn的優(yōu)缺點，經(jīng)過t檢驗篩選后的不同參數(shù)r下的ApEn與SampEn的均值差距率同樣被計算，其值隨r的變化曲線如圖3所示.從圖2、圖3可以看出，r=0.1～0.55，SampEn的均值差距率均高于ApEn的均值差距率，其中經(jīng)過t檢驗篩選后的效果更加明顯，這說明了SampEn方法能更好地區(qū)分HC組與IBS組兩組被試.

圖2 ApEn與SampEn的均值差距率(未經(jīng)過t檢驗篩選)隨r的變化曲線

圖3 ApEn與SampEn的均值差距率(經(jīng)過t檢驗篩選)隨r的變化曲線

2.2 近似熵與樣本熵算法的運行時間分析

本文對被試所有體素的ApEn與SampEn的運行時間進(jìn)行了統(tǒng)計，每個被試所需的計算時間如表3所示.從表3可以看出，相比于ApEn，每一個被試所有體素的SampEn的計算時間平均節(jié)省28秒，這說明SampEn方法的運行效率更高，這一點在大批量fMRI數(shù)據(jù)的分析處理中是很有意義的.

表3 一個被試所有體素的ApEn與SampEn的計算時間

2.3 近似熵與樣本熵方法確定的差異區(qū)域分析

根據(jù)上述計算結(jié)果，可以得到ApEn方法與SampEn方法的最優(yōu)參數(shù)r.根據(jù)Sokunbi等[5]研究中ApEn最優(yōu)參數(shù)r的選取方法，應(yīng)當(dāng)選取較為穩(wěn)定的r值為最優(yōu)參數(shù)，即r的微小變化不會導(dǎo)致全腦ApEn發(fā)生突變，因此，本文ApEn最優(yōu)參數(shù)r取為0.25.為保證ApEn與SampEn方法中參數(shù)r的一致性，SampEn方法的參數(shù)r同樣取為0.25.采用雙樣本t檢驗對平滑后的HC與IBS兩組的ApEn圖及SampEn圖進(jìn)行統(tǒng)計分析，結(jié)果發(fā)現(xiàn)，ApEn差異區(qū)域分布如表4所示，而SampEn差異區(qū)域分布如表5所示.

表4 HC與IBS的ApEn差異區(qū)域分布

表5 HC與IBS的SampEn差異區(qū)域分布

ApEn與SampEn揭示的是大腦功能活動的復(fù)雜度大小，其值越大，表示大腦功能活動越復(fù)雜[9,10].從表4與表5可以看出，IBS患者與HC的ApEn與SampEn均存在顯著性差異，這種差異可能暗示了IBS患者大腦功能活動的異常.本文采用ApEn方法與SampEn方法發(fā)現(xiàn)的IBS患者在補(bǔ)充運動區(qū)、顳區(qū)以及中央旁小葉的異常在以前的IBS或內(nèi)臟疼痛研究[11,12]中均已被證實.

本文采用SampEn方法發(fā)現(xiàn)IBS患者在左側(cè)額上回、右側(cè)額中回及雙側(cè)中央后回的SampEn顯著小于HC，而采用ApEn方法僅發(fā)現(xiàn)IBS患者在右側(cè)中央后回的ApEn顯著小于HC.在以前的研究中，Ma等[13]采用低頻振幅分析方法發(fā)現(xiàn)IBS患者左側(cè)額上回、右側(cè)額中回及雙側(cè)中央后回的功能活動存在異常，與SampEn結(jié)果一致.這暗示了SampEn能夠檢測ApEn無法檢測的腦體素異常.另外，從表4與表5還可以看出，SampEn方法確定的差異簇大小明顯大于ApEn方法.在一定程度上，這說明SampEn方法比ApEn方法更靈敏.

2.4 近似熵與樣本熵算法的理論分析

對比ApEn與SampEn方法在fMRI的實現(xiàn)過程可以發(fā)現(xiàn)，ApEn方法在計算向量之間的距離時包含了自身，而SampEn方法并不包含自身.一方面，由于包含自身與新信息觀點是不相容的，這將導(dǎo)致ApEn方法得到的結(jié)果會存在偏差，因此SampEn方法可以得到更加準(zhǔn)確的結(jié)果，這與本研究中SampEn方法更易區(qū)分HC組與IBS組兩組被試是一致的；另一方面，由于SampEn方法未進(jìn)行自身匹配這一步驟，這解釋了本研究中SampEn方法的運行時間低于ApEn方法的運行時間.綜上所述，SampEn方法可能比ApEn方法更適用于靜息態(tài)fMRI數(shù)據(jù)的分析.

3 結(jié)論

本文針對HC與IBS兩組人的靜息態(tài)fMRI數(shù)據(jù)，采用ApEn方法與SampEn方法進(jìn)行分析，計算不同參數(shù)r下的ApEn與SampEn的均值差距率并進(jìn)行對比.結(jié)果發(fā)現(xiàn)，隨著參數(shù)r的增加，SampEn的均值差距率均比ApEn大.同時，SampEn執(zhí)行單個樣本可節(jié)約28 s的運行時間，這對fMRI大批數(shù)據(jù)的分析處理具有重要的意義.因此，相比于ApEn方法，SampEn方法可能更適用于靜息態(tài)fMRI數(shù)據(jù)的分析.

需要指出本研究存在的不足之處：①本研究采用均值差距率來區(qū)分HC組與IBS組兩組被試的差異性，該指標(biāo)的合理性及有效性仍需更多的fMRI研究來進(jìn)一步驗證；②當(dāng)前已存在多種熵技術(shù)被應(yīng)用于腦研究[14-20]，但本文僅比較了兩種熵在fMRI應(yīng)用中的均值差距離，在將來的研究中，還需要對更多的熵技術(shù)進(jìn)行深入探討.