一種功能核磁共振時間序列匹配方法
【專利摘要】本發明公開了一種功能核磁共振時間序列匹配方法，包括以下步驟：計算系數矩陣β的初始值；計算觀測信號時間序列Y與設計時間序列X間的范式距離Fdist；隨機獲取新的記錄點；計算新的范式距離Fdist；輸出β，完成Y與X的β匹配。本發明在fMRI時間序列匹配中將時域信號經快速離散傅立葉模變換后得到頻域序列，從而完全消除相位信息，以達到消除要匹配的fMRI時間序列間相位差的目的。本發明與當前相位校正的方法相比，更為簡單，消耗的計算量少。本發明對離散傅里葉變換進行加權約束，以降低高頻部分的影響，優先考慮了低頻有效信號，更為明確地確定“最感興趣頻率”的位置，從而降低激活檢驗結果出現“偽陽性”體素的可能性。
【專利說明】一種功能核磁共振時間序列匹配方法
【技術領域】
[0001] 本發明涉及一種數據處理技術，特別是一種功能核磁共振時間序列同頻不同相的 匹配方法。
【背景技術】
[0002] 功能磁共振（fMRI)技術，是無損地研宄人腦行為活動的重要影像技術，具有廣泛 而重要的應用前景。fMRI信號是一種難以直觀解讀的觀測型信號，因此針對fMRI信號的合 理分析具有重要意義，這既是熱點但也是難點問題。
[0003] fMRI技術的通過大腦的血氧水平依賴效應呈現大腦的活動情況，通過該效應，可 以實時地檢測大腦活動的腦血流變化情況。fMRI數據有以下幾個特點：
[0004] 1)延遲特性。一個腦功能刺激，fMRI信號往往并不能做出及時的反饋，而是在幾 秒后達到峰值。
[0005] 2)混淆特性。在以往的信號分析研宄中，fMRI信號通常被看作一種由不同刺激經 過線性疊加而形成的觀測信號。
[0006] 3)采樣點多。fMRI的掃描精度比較高，若以2秒每幀的速度進行掃描，可以獲得 64*64*32分辨率的大腦圖像。因此對維數比較敏感的算法，在應用到fMRI數據中往往會遇 到困難。
[0007] 3)信噪比低。大腦是一個復雜的整體，在進行某些特定認知任務的時候，大腦還 肩負著人體的調節工作，比如呼吸、心跳等一系列維持生命所必須的生理活動。而在進行有 意識的活動的同時，大腦也會出現許多伴隨著諸多因長期的生活訓練而出現的無意識的活 動。
[0008] fMRI實驗數據按照實驗刺激的不同主要分為靜息態數據和功能實驗數據。靜息 態數據無需被掃描者做出任何活動，主要用于觀察人處于無意識狀態下的大腦活動。而功 能數據則是在進行腦科學實驗時，實驗設計者給予受試者外界刺激，誘導受試者有意識地 進行思維活動所產生的fMRI數據，主要用于觀察當人處于有意識的行為時，大腦的活動狀 --τ O
[0009] 目前用于分析fMRI行為實驗數據的常用軟件主要有統計參數映射模型（SPM)、 fMRI組獨立成分工具箱（GIFT)、靜息態數據處理方法（Dparsf)等，這些工具已經用于進行 人腦行為認知的研宄當中，最為廣泛的是SPM(Statistical Parametric Maps)，它可以提 取特定刺激因素引起大腦活動的范圍，即腦激活區。SPM主要利用刺激時間序列這一先驗條 件，提出一種線性疊加模型，即廣義線性模型（GLM)，來分析不同腦區的fMRI時間序列信號 與刺激時間序列的線性疊加關系，從而判定該區域與特定刺激的相關性。SPM框架包括多個 組成部分，包括時間校正、頭動校正、標準化預處理等數據預處理方法，并集成GLM、pair&no pair T檢驗、F檢驗和貝葉斯等多種數據分析方法。Rest、Dparsf、GIFT、MICA等fMRI數據 分析工具包直接集成或加載SPM的部分功能，或要求數據由SPM進行預處理后再進行數據 分析。
[0010] SPM雖然被廣泛應用，但也有部分fMRI信號并沒有得到很好地轉換為腦區激活信 息，因為其依賴的廣義線性模型（SPM-GLM)對刺激時間序列的要求比較苛刻，這大大限制 了 fMRI實驗設計的靈活性，而且該方法將fMRI時間序列視為一個向量，通過評價擬合參考 序列與待測時間序列信號的范數距離作為目標函數，即待測fMRI信號與模型參考時間序 列的匹配程度，來評價優化的目的。但是作為時間序列，相鄰的時間點間聯系是很緊密的， 而fMRI信號作為時間序列，由刺激引發的信號發生延遲的個體差異較大，即使同一個體對 不同的刺激反應延遲也存在差異，因此當腦神經元活動的時間與刺激呈現的時間存在較大 不確定性的差異時，即大腦活動的時間序列本身難以準確描述的時候，傳統的SPM-GLM方 法很難得到理想的結果。
【發明內容】
[0011] 為解決現有技術存在的上述問題，本發明要提出一種簡便的、可以克服延遲的功 能核磁共振時間序列匹配方法。
[0012] 為了實現上述目的，本發明的技術方案如下：一種功能核磁共振時間序列匹配方 法，包括以下步驟：
[0013] A、計算系數矩陣β的初始值
[0014] 將廣義線性模型GLM表述為如下形式：
[0015] Y = Χβ + ε..............(I)
[0016] 其中，Y為觀測信號時間序列，X為設計時間序列，β為系數矩陣，ε為誤差，則系 數矩陣β的初始值的計算公式如下：
[0017] β = (XTX)_1XTy.........(2)
[0018] 式中y是Y的一個分量。
[0019] B、計算觀測信號時間序列Y與設計時間序列X間的范式距離Fdist:
[0020] Fdist (Y, X β ) = norm (Υ-Χ β ) *norm (|fft(Y) |-|fft(X0) | *w).....(3)
[0021] 其中，norm()表示二范數，fft()表示快速傅立葉變換，I · I表示取絕對值，w為 權重向量，其j分量w (j)由下列公式計算：
[0022] w(j) = l/(e_nJ-e0J)...........(4)
[0023] w = w/max (w)............(5)
[0024] 其中n為正向感興趣頻率常數，θ為負向感興趣頻率常數，它們共同通過公式 (5)來確定"最感興趣頻率"的位置。
[0025] C、設置初始值
[0026] 定義初始記錄器Vd= bmin= β，設置ε = 〇. 1，Δ = 1，設置迭代計算器t = 0， sJ= Fdist (y，χ*β )，這里 J· = !。
[0027] D、設置隨機方向正確次數記錄器
[0028] 設置隨機方向正確次數記錄器k = 0。
[0029] E、隨機獲取新的記錄點
[0030] 更新步長Δ = Δ/2,在當前記錄器bold鄰域內，從Δ)分布中隨機獲取 新的記錄點 bnew，使得 bnew滿足（b ^-b-) * ODnew-D + ε >〇。
[0031] F、計算新的范式距離Fdist
[0032] Snew= Fdist (y, x*bnew)............(6)
[0033] 其中x是X的一個分量。
[0034] G、如果sn"〈Sj，則k增加1位。
[0035] H、如果k>v*p，返回步驟D。這里V是采樣數，p是0?0. 5之間的隨機數。
[0036] I、更新 sj+1= s new; β j+1= b min; j 增加 1 位。
[0037] J、轉步驟D，迭代直到j到達迭代最大次數，所述的最大次數為90-100次。
[0038] K、輸出β，完成Y與X的β匹配。
[0039] 與現有技術相比，本發明具有以下有益效果：
[0040] 1、頻域模變換：如公式（3)所示，本發明在fMRI時間序列匹配中將時域信號經快 速離散傅立葉模變換后得到頻域序列，從而完全消除相位信息，以達到消除要匹配的fMRI 時間序列間相位差的目的。本發明與當前相位校正的方法相比，更為簡單，消耗的計算量 少。
[0041] 2、加權約束：考慮離散傅里葉變換，高頻部分的離散傅里葉變換相當于高頻噪聲， 本發明對離散傅里葉變換進行加權約束，如公式（4)-(5)所示，以降低高頻部分的影響，優 先考慮了低頻有效信號，更為明確地確定"最感興趣頻率"的位置，從而降低激活檢驗結果 出現"偽陽性"體素的可能性。
[0042] 3、變步長啟發式尋優：在fMRI中，若將GLM看作一個時間序列的匹配問題，欲尋求 其最優解則無疑是一種連續空間中的問題。如步驟D到步驟I，本發明通過一個隨機給定的 初始點，逐步尋求其鄰域內的更好的點，并利用近似兩個互相穿過對方球心的球的球面所 圍成的一個空間形成的一個二階連續空間作為空間禁忌啟發式準則，提高了搜索效率，當 逐步尋優遇到瓶頸時，逐步減少搜索半徑，提高了搜索精度。
【專利附圖】
【附圖說明】
[0043] 本發明共有附圖4張，其中：
[0044] 圖1是非線性加權向量曲線圖。
[0045] 圖2是本發明的尋優流程圖。
[0046] 圖3是模擬的時間序列匹配圖。
[0047] 圖4是所得的感興區的激活圖。
【具體實施方式】
[0048] 下面結合附圖對本發明進行進一步地描述。
[0049] 圖1是非線性加權向量曲線圖，其中，實線是觀測時間序列變換到頻域后的位置 曲線，虛線是頻域匹配非線性加權權重w。在步驟B中，為fMRI時間序列匹配提供頻域匹配 權重w，也就是圖1中虛線所示，提高信噪比，降低激活檢驗結果出現"偽陽性"體素的可能 性。
[0050] 圖2是本發明的尋優流程圖，該流程圖是在時間序列頻域對齊之后通過迭代方式 降低匹配誤差的過程，如步驟D到J所描述，該尋優方式是在二階禁忌空間中進行隨機鄰域 優化，使得本發明中的匹配可以穩定收斂到全局最優解，同時避免陷入局部最優解，而且如 步驟E所描述，尋優步長△可變，提高了本發明的精度，使得最終得到的β解可以穩定在 精度范圍內，并且根據Abel收斂準則，該收斂于全局最優解。
[0051] 圖3是模擬的時間序列匹配圖，其中實線是設計時間序列X，帶小方塊標記的實線 是由本發明進行匹配所得到的fMRI信號時間序列，顯示了本發明的fMRI時間序列匹配效 果。
[0052] 圖4是所得的感興區的激活圖，該圖是經過本發明中步驟A到K完成記憶想象實 驗觀測時間序列與設計時間序列匹配后對腦激活的判別結果，從圖4可以看出激活體素準 確。
[0053] 表1是激活體素數目統計結果，從比較結果看，經過本發明的匹配方法檢測激活 體素，有效降低"偽陽性"。
[0054] 表1 :激活體素數目統計結果表
[0055]
【權利要求】
1. 一種功能核磁共振時間序列匹配方法，其特征在于：包括以下步驟： A、 計算系數矩陣0的初始值 將廣義線性模型GLM表述為如下形式： Y = X0 + e ............... (1) 其中，Y為觀測信號時間序列，X為設計時間序列，0為系數矩陣，e為誤差，則系數矩 陣0的初始值的計算公式如下： 0 = (XTX)^XTy......... (2) 式中y是Y的一個分量； B、 計算觀測信號時間序列Y與設計時間序列X間的范式距離Fdist: Fdist (Y, X 0 ) = norm(Y-X 0 ) *norm( | fft (Y) | -1 fft (X 0 ) | *w) --?. . (3) 其中，norm()表示二范數，fft()表示快速傅立葉變換，| ? |表示取絕對值，w為權重 向量，其j分量w(j)由下列公式計算： w(j) = l/(e-nJ-e0J).........-- (4) w = w/max (w)............ (5) 其中n為正向感興趣頻率常數，0為負向感興趣頻率常數，它們共同通過公式（5)來 確定"最感興趣頻率"的位置； c、設置初始值 定義初始記錄器13。1(1= bmin= |3，設置e = 〇. 1，A = 1，設置迭代計算器t = 0，s」= Fdist (y，x* 0 )，這里 j = 1 ; D、 設置隨機方向正確次數記錄器 設置隨機方向正確次數記錄器k = 0 ; E、 隨機獲取新的記錄點 更新步長A = A/2,在當前記錄器bold鄰域內，從A)分布中隨機獲取新的 記錄點 bnew，使得 bnew滿足（b min-bQld) * (bnew-bQld) + e >〇 ; F、 計算新的范式距離Fdist Snew= F dist(y? X*bnew)............ (6) 其中X是X的一個分量； G、 如果snOT〈Sj，則k增加1位； H、 如果k>v*p，返回步驟D ;這里v是采樣數，p是0?0. 5之間的隨機數； I、 更新 sj+1= s new; |3 j+1= b min; j 增加 1 位； J、 轉步驟D，迭代直到j到達迭代最大次數，所述的最大次數為90-100次； K、 輸出0，完成Y與X的0匹配。
【文檔編號】A61B5/055GK104434109SQ201410805041
【公開日】2015年3月25日 申請日期:2014年12月19日 優先權日:2014年12月19日
【發明者】劉洪波, 陳亮, 張維石, 馮士剛 申請人:大連海事大學