楊增宇，蔣文濤，陳宇

（四川大學工程力學系，成都610065）

1 引 言

近年來，功能磁共振成像（fMRI）已經(jīng)成為研究大腦動作方式的首要選擇方法［1-3］。其原理主要是，BOLD-FMRI圖像以體素為單位把神經(jīng)元活動引起的生理活動記錄下來，針對每個體素，我們觀測得到的FMRI數(shù)據(jù)相當于一個時間序列，如果將神經(jīng)元的活動引起的生理活動的變化統(tǒng)稱為神經(jīng)活動的血液動力學響應，那么可以采用如下的方式來描述：觀測到的體素時間序列＝血液動力學響應＋漂移＋噪聲［4］。因此，要探究fMRI數(shù)據(jù)與神經(jīng)活動的關系，研究血液動力學響應是關鍵。如果用X（t）表示血流動力學響應，u（t）表示神經(jīng)活動，h（t）表示相應的血液動力學響應函數(shù)，那么在線性時不變系統(tǒng)模型［5-6］中 X（t），u（t），h（t）之間的關系可以表示為：X（t）＝u（t）×h（t），即血流動力學響應等于神經(jīng)活動與血液動力學響應函數(shù)的卷積。因此一旦對血液動力學響應函數(shù)有了相當?shù)恼J識就可以預測由某種神經(jīng)活動模式所導致的血流動力學響應進而得到體素時間序列，實現(xiàn)fMRI成像功能。

2 標準血液動力學響應模型

由于血液動力學響應函數(shù)直接影響觀測的體素時間序列，所以對于血液動力學響應函數(shù)的研究具有重要意義。目前三參數(shù)gamma函數(shù)形式是fMRI數(shù)據(jù)處理的研究中應用最為廣泛的血液動力學響應函數(shù)之一。1994年Friston等［7］提出可以用gamma函數(shù)的形式來表征血液動力學響應函數(shù)，在此基礎上Cohen［8］于1996年提出了三參數(shù)形式的血液動力學響應函數(shù)模型：

眾多研究［9］已經(jīng)表明上述標準的血液動力學響應函數(shù)模型在描述和預測血液動力學響應時具有一定的合理性和有效性，然而由于在實驗數(shù)據(jù)擬合時未考慮參數(shù)間的耦合作用，導致得到的函數(shù)模型在求解精度和擬合效果上都并非最優(yōu)，迄今國內外針對血液動力學響應函數(shù)的研究幾乎沒有對該三參數(shù)模型進行相應的優(yōu)化處理工作。本研究將以此模型為基本框架，通過擬合實驗數(shù)據(jù)值對其進行優(yōu)化。

3 參數(shù)敏感性分析

由3參數(shù)控制的血液動力學響應模型，由于各參數(shù)對最終信號響應值的影響程度是不同的：影響大的參數(shù)即使有微小攝動，都會引起流動應力的較大波動，甚至吞沒了其他相對不敏感參數(shù)的作用。因此，有必要對由實驗數(shù)據(jù)確定的血液動力學響應模型的參數(shù)進行敏感度分析。

3.1 參數(shù)敏感度分析方法

拉丁超立方抽樣（LHS）是一種均勻空間抽樣方法，最早于1979年由simpson等提出，具有抽樣次數(shù)較少、可有效避免重復抽樣等特點，適用于對復雜的多維參數(shù)空間抽樣［10］。抽樣過程如下：

（1）設變量X的取值范圍為［a，b］，所需抽取的樣本個數(shù)為n。

（2）將X的取值范圍［a，b］等概率的分為 n個區(qū)間，并在每個區(qū)間中隨機的抽取一個數(shù)作為Xi的一個樣本值，將所取的值進行隨機排列，由此可得到關于Xi的樣本容量為n的樣本序列：｛x1，x2，x3……xn｝。

（3）對于m維參數(shù)空間（即含有m個變量）的抽樣，則是對 m個變量（X1，X2，X3……Xm）均進行步驟1和步驟2，得到m個樣本序列。

（4）對m個樣本序列各取一個值即可獲得一個m維的樣本參數(shù)集，共可得到n個m維的樣本參數(shù)集，見表1，此即完成了拉丁超立方抽樣的全過程。

表1 多參數(shù)LHS抽樣樣本參數(shù)集Table 1 Sample parameter set of Multi parameter LHS sampling

Spearman秩相關系數(shù)是一個非參數(shù)性質（與分布無關）的秩統(tǒng)計參數(shù)［11］，由 Spearman在1904年提出，用于度量兩個隨機變量之間的相關方向與密切程度，在統(tǒng)計學中，Spearman秩相關系數(shù)一般由希臘字母ρs表示。設變量X的樣本集為｛x1，x2，x3……xn｝，將其按照從大到小的順序進行排列，得到一個有序序列｛x1′，x2′，x3′……xn′｝。定義原始 xi在新的有序序列｛x1′，x2′，x3′……xn′｝中所占的位置Ki為 xi在樣本集｛x1，x2，x3……xn｝中的秩。設與 X對應的變量Y的樣本參數(shù)集為｛y1，y2，y3……yn｝。設 xk在｛x1，x2，x3……xn｝中的秩為 αk，yk在｛y1，y2，y3……yn｝中的秩為 βk，則相應的 spearman秩相關系數(shù)的計算公式為：

Spearman秩相關系數(shù)即反映了該參數(shù)的敏感度［8］。

3.2 參數(shù)敏感度分析結果

對模型參數(shù)（k、n、m）進行超立方抽樣，形成 p個3維的樣本參數(shù)集，見表2，每個參數(shù)集對應一個具體血液動力學響應函數(shù)模型。

表2 血液動力學響應模型的LHS抽樣樣本參數(shù)集Table 2 LHS sampling parameter set of hemodynamic response model

將所得到的各個樣本參數(shù)集的參數(shù)帶入模型中，依據(jù)文獻［8］中實驗的時間節(jié)點值計算出相應的信號強度h*，而文獻［8］中實驗測得的信號強度為h，根據(jù)最小二乘擬合的原理，計算目標函數(shù)值W，依據(jù)spearman秩的定義和式（1），利用所得的 W值分別和相應的 k、n、m值計算spearman秩相關系數(shù)，從而得到參數(shù)敏感度，見表3。

表3 模型參數(shù)敏感度統(tǒng)計表Table 3 Statistics of model parameter sensitivity

由表3中數(shù)據(jù)可知：（1）當各參數(shù)的取值范圍一定時，響應函數(shù)的參數(shù)敏感性估計受抽樣次數(shù)的影響，并且抽樣次數(shù)越大，各參數(shù)的敏感度越趨于穩(wěn)定。所以為了取得較真實、較準確的模型參數(shù)敏感度，應將抽樣次數(shù)取得足夠大，以保證所得到的敏感性分析結果有效性及合理性。當抽樣次數(shù)由500次增加至1 000次時參數(shù)敏感度波動很小，故可認為500次抽樣結果已經(jīng)足夠穩(wěn)定，能夠代表樣本特征，此時的數(shù)據(jù)具有分析和應用價值。（2）當保持參數(shù)k、m取值范圍以及抽樣次數(shù)一定（均為500次），只改變參數(shù)n的取值范圍時，各參數(shù)敏感度均存在較大的波動，故各參數(shù)敏感度不僅與其自身的取值范圍有關，還與其他參數(shù)的取值范圍相關，即模型參數(shù)間的耦合作用對分析結果的影響很大。（3）三個參數(shù)k、m、n中，m和 n的參數(shù)敏感度相當，而 k的參數(shù)敏感度較小，與m、n相比有很大差異。即m、n對目標函數(shù)的影響較大，為了使各個參數(shù)都有相當?shù)膮?shù)敏感度，應該使各個參數(shù)對目標函數(shù)的spearman秩相關系數(shù)都接近0.333，然而這樣又會將k值的取值區(qū)間設置得過小，不利于后續(xù)的優(yōu)化搜索工作，故我們取表中k、m、n參數(shù)敏感度最接近的情況，即取初始取值區(qū)間為：k：［0，10］m：［0，40］n：［-20，20］

4 優(yōu)化與分析

4.1 模型優(yōu)化算法—遺傳算法

遺傳算法是一種模仿生物的遺傳進化原理，通過選擇、交叉與變異等操作機制，使種群中個體的適應性不斷提高，從而找到問題的滿意解或最優(yōu)解的仿生全局優(yōu)化算法［12］。

遺傳算法基本流程框見圖1。在參數(shù)的初始取值區(qū)間范圍內隨機取值，作為一個個體，將所有初始個體組成的群體稱為初始種群，對初始種群中每個個體進行二進制編碼，所得到的二進制字符串通常稱為該個體的“染色體”；設置合理的適應度函數(shù)，計算初始種群的適應度值，并進行是否終止進化的判斷；對于種群需要多次進化才能滿足求解要求的問題，對上一代種群進行選擇（輪盤賭方法選擇新一代進化種群）、交叉（兩兩個體之間按照交叉概率交換他們的部分染色體片段）、變異（每個個體按照變異概率改變某一個或某些基因座上的基因值為其它的等位基因）等運算得到新一代種群，并判斷是否滿足終止進化條件。

4.2 模型優(yōu)化分析

根據(jù)參數(shù)敏感性分析結果，對模型中的3個參數(shù)給定初始取值區(qū)間：k：［0，10］，m：［0，40］，n：［-20，20］，利用matlab平臺編寫遺傳算法程序，以文獻［8］中實驗值為基礎，分別取遺傳代數(shù)為50、100、500、1 000對模型參數(shù)進行最優(yōu)值求解，最終，求得使目標函數(shù)W（x）最小時所對應的參數(shù)即響應模型：

圖1 遺傳算法基本流程框圖Fig 1 Basic flow diagram of genetic algorithm

50代：［k，m，n］＝［21.654，1.878，-0.891］h（t）＝21.654t1.878e-0.891t

100代：［k，m，n］＝［17.017，7.393，-2.358］h（t）＝17.017t7.393e-2.358t

500代：［k，m，n］＝［2.095，6.699，-1.538］h（t）＝2.095t6.699e-1.538t

1000代：［k，m，n］＝［0.518，7.780，-1.601］h（t）＝0.518t7.780e-1.601t

2000代：［k，m，n］＝［0.457，8.460，-1.804］h（t）＝0.457t8.460e-1.804t

將所得模型與文獻［8］中實驗值進行擬合效果比較：

圖2 不同遺傳代數(shù)擬合效果比較Fig 2 Comparison of fitting effect inthe different genetic algebra

由圖2可知隨著遺傳代數(shù)的增加，通過優(yōu)化搜索所得的模型的擬合效果越來越好；算法顯示出較快的收斂性，當遺傳代數(shù)為100代時所得模型已經(jīng)初步具有實驗值曲線的輪廓；當遺傳代數(shù)多于1 000代時，代與代之間差距較小，已經(jīng)接近最優(yōu)的模型曲線。

4.3 模型優(yōu)化結果

取遺傳代數(shù)為2500代（可以認為此時目標函數(shù)已收斂到最小值），進行遺傳算法最優(yōu)搜索，得到優(yōu)化的模型為：

而根據(jù)文獻［8］中實驗值，由Cohen利用最小二乘法確定的三參數(shù)形式的血液動力學響應函數(shù)模型為：

圖3 Cohen模型與本研究模型的比較Fig 3 Comparison of the Cohen model with the model presented in this study

由圖3可知Cohen模型與本研究所建立的模型都能夠描述相應的血流動力學響應，但是本研究所建立的模型能夠更加貼合實驗曲線，具有更高的求解精度與更好的擬合效果。

5 經(jīng)典算例分析

為了說明優(yōu)化的響應模型能夠較好的描述與預測由神經(jīng)活動（u）所引起的血液動力學響應（X），本研究以一個OFF-ON任務為例，任務數(shù)據(jù)取自文獻［13］Savoy及其學生做的視覺閃光實驗，實驗采用方波刺激，有交替的45 s的刺激和休息時間，實驗記錄下fMRI信號隨時間的變化過程。利用Cohen模型以及所得到的優(yōu)化后的響應模型式（4）來進行比較與驗證，可以得到圖4的血液動力學響應，圖中灰色區(qū)域表示刺激期，白色區(qū)域表示休息期，連續(xù)曲線為神經(jīng)活動（此處為δ（t）階躍函數(shù)）和響應函數(shù)的卷積的形成的響應預測曲線?？梢钥吹?，（1）響應信號相對任務的開始時間有一定的延遲，在任務結束后也需要一定的時間回到基線水平，這與通常觀測到的fMRI信號特點是吻合的。（2）與傳統(tǒng)的Cohen模型相比較，本研究所得優(yōu)化模型能夠更好描述和預測由神經(jīng)活動所引起的血液動力學響應，所得結果更加接近真實值。

圖4 OFF-ON任務響應圖Fig 4 OFF-ON task response diagram

6 結論

血液動力學響應函數(shù)模型參數(shù)的敏感度受其他參數(shù)的取值范圍的影響，故不能單獨研究某一材料參數(shù)的敏感度，表明本研究的LHS和spearman秩相關結合的整體分析方法對于該問題是適合的。本研究采用參數(shù)敏感度估計搜索區(qū)域和遺傳算法結合的方法能夠很好的實現(xiàn)血液動力學響應函數(shù)模型的參數(shù)識別。所用于最優(yōu)化處理的遺傳算法具有較好的收斂性與普適性。建立的優(yōu)化模型能夠很好的描述和預測由神經(jīng)活動引起的血流動力學響應，所得到的預測曲線與通常觀測到的fMRI數(shù)據(jù)特點吻合。此外，本研究建立的優(yōu)化模型在求解精度和擬合效果上較傳統(tǒng)的Cohen模型有一定的提高。