劉 婷,戴亞康,楊瑩雪,王玉平

(1.中國科學院蘇州生物醫(yī)學工程技術(shù)研究所,江蘇蘇州 215163;2.中國科學院長春光學精密機械與物理研究所,吉林長春 130033;3.中國科學院大學,北京 100049;4.首都醫(yī)科大學宣武醫(yī)院神經(jīng)內(nèi)科,北京 100053;5.腦功能疾病調(diào)控治療北京市重點實驗室,北京 100053 )

基于時域平滑約束的腦磁時序信號逆問題求解方法

劉 婷1,2,3,戴亞康1,楊瑩雪4,5,王玉平4,5

(1.中國科學院蘇州生物醫(yī)學工程技術(shù)研究所,江蘇蘇州 215163;2.中國科學院長春光學精密機械與物理研究所,吉林長春 130033;3.中國科學院大學,北京 100049;4.首都醫(yī)科大學宣武醫(yī)院神經(jīng)內(nèi)科,北京 100053;5.腦功能疾病調(diào)控治療北京市重點實驗室,北京 100053 )

由腦磁時序信號重建腦內(nèi)時序神經(jīng)信號時,除了要保證重建信號位置和強度的準確性,還要避免重建源信號在時域上瞬變.針對這一問題,提出了一種基于時域平滑約束的腦磁時序信號逆問題求解方法.該方法不同于傳統(tǒng)最小范數(shù)估計算法(Minimum Norm Estimate,MNE),通過引入時域平滑正則算子構(gòu)造雙參數(shù)混合正則化,根據(jù)廣義交叉驗證(Generalized Cross-Validation,GCV)原則選取雙正則化參數(shù)后,根據(jù)單正則項的解在源信號中的權(quán)重將其進行線性組合估算出源信號.仿真數(shù)據(jù)實驗表明,本文方法比傳統(tǒng)MNE方法的總體均方誤差小,且各時刻均方誤差基本穩(wěn)定在同一水平;同時本文方法重建的源信號與仿真源信號變化趨勢基本一致.真實數(shù)據(jù)實驗發(fā)現(xiàn),本文方法重建結(jié)果的曲率變化率為0.0640,而傳統(tǒng)MNE方法重建結(jié)果的曲率變化率為0.1646.實驗結(jié)果證明本文方法能重建出空域準確且時域平滑的腦內(nèi)神經(jīng)信號.

腦磁時序信號;逆問題;雙參數(shù)混合正則化;時域平滑

1 引言

人類的腦部活動與大腦細胞活動息息相關(guān),其中最主要的腦細胞——神經(jīng)元——是大腦神經(jīng)傳遞的主體,它通過產(chǎn)生動作電位來傳遞神經(jīng)信號.腦物理學認為,如果把人類大腦看作是電磁系統(tǒng)時,它遵守物理學中的電磁規(guī)律,由受到刺激或心理活動激活的神經(jīng)元充當激勵腦電磁場的源[1].腦磁圖(MagnetoEncephaloGraphy,MEG)是測量腦神經(jīng)信號的非侵入性腦功能檢測技術(shù),它利用超導線圈在腦外測量神經(jīng)元產(chǎn)生的微弱磁場,通過分析MEG信號反推大腦內(nèi)部神經(jīng)元活動,此即腦磁逆問題.腦磁逆問題有兩個難點:一是解的非唯一性,二是解的不穩(wěn)定性.為了能夠得到穩(wěn)定、合理的解,必須通過先驗知識在逆問題求解過程中對解空間加以約束,即引入正則化技術(shù)[2,3].從上個世紀90年代開始科學家們提出了諸多腦內(nèi)源定位技術(shù),其中開展最早且被普遍采用的方法是最小范數(shù)估計算法[4](Minimum Norm Estimate,MNE).此算法的本質(zhì)是尋找最小能量解,通過采用Tikhonov正則化[3]平衡數(shù)學模型誤差與源能量,以全局能量最小的源來推算腦內(nèi)源信號的位置和強度.

近年來隨著研究的進展,人們開始探索諸如注意機制、沖突處理等的神經(jīng)傳導機制.科學家們不再滿足于只獲知單一時刻腦內(nèi)源信號的位置、強度,還希望獲知腦內(nèi)神經(jīng)信號的定向傳導過程,也就是腦磁源時序信號.腦神經(jīng)學認為,相鄰腦結(jié)構(gòu)之間的興奮傳導間隔幾個到幾十個毫秒之間[5],而不會有相鄰時刻間跳變的現(xiàn)象.由于沒有施加時域約束,傳統(tǒng)MNE方法重建出來的腦磁源信號往往偏離原始時序源信號的變化趨勢,在時域上有明顯的跳變現(xiàn)象,這與神經(jīng)信號的定向傳導機制相違背[6].為此,本文提出基于時域平滑約束的MEG時序信號逆問題求解方法.該方法從傳統(tǒng)MNE方法出發(fā),在Tikhonov正則化中引入時域平滑約束項,構(gòu)造雙參數(shù)混合正則化,根據(jù)廣義交叉驗證(Generalized Cross-Validation,GCV)準則選取兩個合適正則化參數(shù)后,通過計算單正則項的解在源信號中占的權(quán)重然后進行線性組合估算出源信號.仿真實驗表明,本文提出的方法不僅可以準確重建腦磁源,而且重建的時序源信號能更好地還原真實源信號的變化趨勢,大大改善傳統(tǒng)MNE方法重建結(jié)果在時域上振蕩的現(xiàn)象.該方法將在正文第二部分詳細描述,第三部分是仿真、真實數(shù)據(jù)實驗和實驗結(jié)果,最后是本文的討論部分.

2 理論和方法

假設(shè)腦外有m個通道的MEG信號,腦內(nèi)有n個均勻分布的源信號,那么在i時刻腦內(nèi)源信號與MEG信號的關(guān)系可以用以下離散化線性模型[7]表示:

bi=Axi+ei

(1)

其中bi為i時刻大小為m×1的MEG測量信號;xi為i時刻腦內(nèi)源信號,大小為n×1;ei是i時刻和bi同維度的噪聲信號;A為轉(zhuǎn)換矩陣,代表腦內(nèi)源信號與MEG測量信號的映射關(guān)系,大小為m×n.矩陣A可以通過邊界元方法、有限元方法等結(jié)合頭模型求解.腦磁逆問題的求解面臨兩個難點:一是解的非唯一性,由于MEG信號通道數(shù)m遠小于腦皮層網(wǎng)格數(shù)n,所以式(1)是高度欠定方程,有無數(shù)個解;二是解的不穩(wěn)定性,即病態(tài)性,由于矩陣A的條件數(shù),即最大特征值與最小特征值之比很大,MEG測量信號中很小的噪聲都將對解產(chǎn)生很大的擾動.基于以上問題,對式(1)求解xi轉(zhuǎn)化為求解最小二次泛函的問題,并且引入Tikhonov正則化技術(shù)來使病態(tài)問題適定化.具體地,在第i時刻,腦磁逆問題求解轉(zhuǎn)化為求解以下最小值問題:

(2)

等式右邊第一項表示測量數(shù)據(jù)和估計數(shù)據(jù)的擬合,第二項為正則項,表示解的先驗信息,其中R為約束解空間的正則算子,λ為正則化參數(shù),調(diào)節(jié)擬合項和正則項在兩項之間達到平衡.當R為單位陣時,式(2)為Tikhonov零階正則化,約束項使解具有全局最小能量;當R為一階或者二階微分矩陣時,式(2)為廣義Tikhonov正則化,約束項使解具有光滑的曲面梯度或曲率.由于矩陣A為行滿秩,計算其Moore-Penrose右逆矩陣,式(2)對應的解的形式為:

(3)

如前所述,式(2)所示目標函數(shù)重建出來的源信號各個時刻之間是相獨立的(可以通過第二部分仿真實驗的結(jié)果看到).為保證時域平滑性,我們對式(2)增加時域平滑約束項,以構(gòu)造雙參數(shù)混合正則化來重建MEG時序源信號,下面介紹具體方法.

2.1 目標函數(shù)構(gòu)造

首先將式(1)轉(zhuǎn)化成時序信號形式.對于時長為k的MEG測量信號,對應的離散化線性模型為:

b=Ax+e

(4)

其中b為MEG測量信號,大小為m×k;A是m×n維轉(zhuǎn)換矩陣,x為n×k維時序源信號矩陣;e為m×k維噪聲信號矩陣.

根據(jù)MNE算法的思想,新的目標函數(shù)既要滿足重建的源信號在整個k時段所有解中能量最小,又要滿足在相鄰時刻間是平滑的[8],因此引入時域平滑約束項,構(gòu)造以下目標函數(shù):

(5)

λ1和λ2都是正則化參數(shù),式(6)右邊第二項作為能量約束項,將逆問題的解約束為在時段k內(nèi)全局能量最小的解;第三項作為時域平滑約束項,使解在相鄰時刻間變化率最小.本文方法具有獨立地對能量和時域平滑約束的性質(zhì):當λ1趨于零時,式(6)主要是時域平滑約束發(fā)揮作用;當λ2趨于零時,式(6)主要是能量約束發(fā)揮作用.λ2為零時本文方法等同于傳統(tǒng)MNE方法,也就是說傳統(tǒng)MNE方法是本文方法的一個特例.在求解時,通過適當?shù)卣{(diào)整正則化參數(shù)λ1和λ2來達到能量約束和時域平滑約束之間的平衡,進而重建出能量小且時域平滑的信號.

Aliev B[9,10]引入類似本文的雙參數(shù)正則化求解線性不適定算子方程時,從數(shù)學的角度證明了普適的雙參數(shù)正則化解的唯一性、穩(wěn)定性和收斂性.王文娟等[11]采用雙參數(shù)正則化方法研究電導率反演成像中發(fā)現(xiàn)雙參數(shù)正則化方法增強了反演的穩(wěn)定性.另外,增加時域平滑約束后,重建信號對正則化參數(shù)的敏感度降低.Brooks等[8]在心電逆問題求解中通過實驗發(fā)現(xiàn),在正則化參數(shù)變化幅度相同的情況下,相比于僅有能量約束,增加時域平滑約束的雙參數(shù)正則化方法重建的信號更穩(wěn)定.因此本文方法具有良好的魯棒性,下文介紹正則化參數(shù)選擇策略和具體求解方法.

2.2 求解方法

首先根據(jù)Kronecker積的定義[12]將式(6)轉(zhuǎn)化成如下形式:

(6)

(7)

3 實驗

3.1 仿真實驗

在腦皮層選取兩個活化位置,坐標分別為(-39.4982,-36.6656,56.8917)和(36.0071,-18.8000,58.9000),兩個位置分別對應左腦和右腦感覺區(qū).6ms時達到能量峰值的源信號被放置在(-39.4982,-36.6656,56.8917)處,19ms時達到能量峰值的源信號則被放置在(36.0071,-18.8000,58.9000)處,圖2所示為未添加噪聲信號時,在第6ms和第19ms時仿真信號在腦皮層和測量空間的成像圖.

本文在用于腦電/腦磁信號分析的開源軟件eConnectome[14,15]平臺上完成了上述仿真數(shù)據(jù)的設(shè)計,并在此基礎(chǔ)上對本文提出的方法進行了實驗驗證.具體地,我們將MEG仿真數(shù)據(jù)導入eConnectome后執(zhí)行數(shù)據(jù)預處理(preprocessing),包括baseline correction (以1～4ms為基準線)和filtering(50Hz陷波濾波器),采用真實幾何頭模型和邊界元方法求解正問題獲取轉(zhuǎn)換矩陣A,然后分別用式(2)傳統(tǒng)MNE方法和式(6)時域平滑約束方法對經(jīng)過預處理的數(shù)據(jù)進行腦磁源重建.

仿真實驗結(jié)果考察兩個方面:一是考察數(shù)據(jù)精確度參數(shù)均方誤差,二是考察兩個活化位置的估算信號與原始模擬信號的吻合情況.

圖4分別展示了腦皮層上兩個活化位置(-39.4982,-36.6656,56.8917)和(36.0071,-18.8000,58.9000)采用傳統(tǒng)MNE方法的估算信號和仿真信號之間的吻合情況.發(fā)現(xiàn)各個時刻間獨立求逆使得解在時域上不規(guī)則振蕩,且某些時刻與真實值相去甚遠.圖5顯示引入雙參數(shù)正則化增加時域平滑約束項后,估算信號基本復原了仿真信號變化趨勢,而且分別在6ms和19ms處具有能量峰值.需要注意的是,估算信號的幅度小于真實信號,是因為式(9)中第二項是能量約束項,也就是說所求的估算信號是所有解中能量最小的解,這是重建算法本身決定的,MNE算法也存在同樣的現(xiàn)象.

3.2 真實數(shù)據(jù)實驗

本實驗數(shù)據(jù)來自首都醫(yī)科大學宣武醫(yī)院神經(jīng)內(nèi)科對焦慮患者的認知實驗數(shù)據(jù)集.研究發(fā)現(xiàn),大腦對沖突信息進行加工時會在刺激出現(xiàn)后的270ms左右誘發(fā)負性相關(guān)電位,即認知電位沖突性負波N270[16].目前的研究認為由N270反映的認知沖突處理系統(tǒng)可能分散在大腦多個不同區(qū)域,但額內(nèi)側(cè)扣帶回可能是該系統(tǒng)的重要組成部分[17].本實驗通過分析一例焦慮患者的認知實驗數(shù)據(jù)來驗證時域平滑約束算法的有效性.實驗為比較圖形的顏色和形狀,每一刺激對中隨機呈現(xiàn)不同的顏色或形狀讓受試者判斷,兩個刺激各自持續(xù)500ms,兩個刺激之間間隔200ms,每個刺激對間隔2s.

實驗通過306導型號為ElektaNeuromag腦磁圖儀采集MEG數(shù)據(jù),分辨率為1000Hz.數(shù)據(jù)經(jīng)過去眼電、濾波和基線校準等預處理步驟后,按照刺激對數(shù)疊加平均成時長為2000ms的MEG數(shù)據(jù).為了定位N270認知沖突處理系統(tǒng),選取1106ms～1285ms(對應N270時序段)為分析時程(epoch),如圖6所示.以額內(nèi)側(cè)扣帶回為興趣區(qū)域,隨機選擇該區(qū)域內(nèi)的某一位置,分別用式(3)傳統(tǒng)MNE方法和式(11)時域平滑約束方法進行源時序信號重建.圖7(a)所示為左前額內(nèi)側(cè)中心坐標為(-3.8603,59.3580,17.9440)處的重建結(jié)果.

參考圖6(b)中的全局能量譜(GlobalFieldPower,GFP)曲線,本文提出的時域平滑約束方法較好地復原了源時序信號的形狀,也成功地定位出兩個能量峰值;而傳統(tǒng)MNE方法在兩個能量峰值處出現(xiàn)了不同程度的抖動,時域平滑性低于本文提出的時域平滑約束方法.我們還通過圖7(b)所示的各時刻曲率變化率(curvaturevariability)定量比較了兩個重建信號的平滑程度,傳統(tǒng)MNE方法的總曲率變化率為0.1646,時域平滑約束方法的總曲率變化率為0.0640,也就是說,對于這個實驗,在重建平滑能力上,時域平滑約束方法要比傳統(tǒng)MNE方法強2.5倍以上.

由上述實驗結(jié)果可知,本文提出的基于時域平滑約束的雙參數(shù)MEG時序信號逆問題求解方法要優(yōu)于傳統(tǒng)MNE方法.

4 討論

大腦對外部信息的反應是復雜的神經(jīng)動力學過程,MEG逆問題的求解從MNE算法的提出雖然已有了很多發(fā)展,但仍有許多問題值得探索.本文從MNE算法出發(fā),引入雙參數(shù)混合正則化,通過時域平滑算子約束解空間,使得估算的時序源信號更符合神經(jīng)信號定向傳導的性質(zhì).這也為因果性腦網(wǎng)絡[18]研究提供了更有效的分析工具.由于MNE算法的缺陷,時域約束估算值不可避免的有“模糊效應”[6],下一步的工作將在此基礎(chǔ)上嘗試減小這一效應,使得腦磁源重建結(jié)果更加精確.

[1]吳殿鴻,郭立文,等.腦物理學[M].哈爾濱:哈爾濱工業(yè)大學出版社,1995.39-41.

[2]肖庭延,于慎根,等.反問題的數(shù)值解法[M].北京:科學出版社,2003.18-30.

[3]Tikhonov A,Arsenin V.Solutions of Ill-posed Problems [M].Washington DC:Winston,1977.113-135.

[5]Fredric MH,Ivica K.神經(jīng)計算原理[M].葉世偉,王海娟,譯.北京:機械工業(yè)出版社,2007.7-11.

[6]Tian S T,Huang J Z,Shen H,Li Z.A two-way regularization method for MEG source reconstruction[J].The Annals of Applied Statistics,2012,6(3):1021-1046.

[7]Bolstad A,Veen B V,Novak R.Space-time event sparse penalization for magneto-/electroencephalography[J].NeuroImage,2009,46(4):1066-1081.

[8]Brooks D H,Ahmad G F,MacLeod R S,et al.Inverse electrocardiography by simultaneous imposition of multiple constraints[J].IEEE Trans.Biomed,1999,46(1):3-81.

[9]Alive B.Two-parameter regularization method for finding L-pseudo-solutions[J].Vestnik Moskov Univ Vychisl Mat Kibernet,1986,15(2):45-50.

[10]Alive B.Modification of the generalized discrepancy principle for L-pseudo-solutions in the degenerate case[J].Vestnik Moskov Univ Vychisl Mat Kibernet,1991,20(1):28-33.

[11]王文娟，Chris Farmer,等.雙參數(shù)混合正則化方法及在電導率反演成像中的應用[J].地球物理學報，2011,54(8)：2154-2159. WANG Wen-juan， Farmer C,et al.A dual-parameter regularization method for electrical conductivity imaging[J].Chinese Journal of Geophysics,2011,54(8):2154-2159.(in Chinese)

[12]張賢達.矩陣分析與應用[M].北京：清華大學出版社,2004.107-117. ZHANG Xian-da.Matrix Analysis and Applications[M].Beijing:Tsinghua University Press,2004.107-117.(in Chinese)

[13]朱南海,趙曉華.基于遺傳算法的Tikhonov正則參數(shù)優(yōu)化計算[J].工程力學,2009,26(5):25-30. ZHU Nan-hai,ZHAO Xiao-hua.Optimal calculation of Tikhonov regularization parameter based on genetic algorithm[J].Engineering Mechanics,2009,26(5):25-30.(in Chinese)

[14]He B,Dai Y K,et al.eConnectome:A MATLAB toolbox for mapping and imaging of brain functional connectivity[J].Journal of Neuroscience Methods,2011,195(2):261-269.

[15]Dai Y K,Zhang W B,et al.Sourceconnectivity analysis from MEG and its application to epilepsy source localization[J].Brain Topography,2012,25(2):157-166.

[16]歐陽取平,王玉平.工作記憶對沖突性負波N270的影響[J].臨床神經(jīng)電生理學雜志,2008,17(6):323-327. OUYANG Qu-ping,WANG Yu-ping.The effects of working memory on the event-related potential N270[J].Journal of Clinical Electroneurophysiology,2008,17(6):323-327.(in Chinese)

[17]王玉平.事件相關(guān)電位N270的特性及本質(zhì)[J].臨床神經(jīng)電生理學雜志,2002,11(4):247-248. WANG Yu-ping.The character and nature of the event-ralated potential N270[J].Journal of Clinical Electroneurophysiology,2002,11(4):247-248.(in Chinese)

[18]孫俊峰,洪祥飛,童善保.復雜腦網(wǎng)絡研究進展—結(jié)構(gòu)、功能、計算與應用[J].復雜系統(tǒng)與復雜性科學,2010,7( 4):74-90. Sun J F,Hong X F,Tong S B.A survey of complex brainnetworks:structure function computation and applications[J].Complex Systems and Complexity Science,2010,7(4):74-90.(in Chinese)

劉 婷 女,1986年出生于江蘇贛榆,現(xiàn)為中科院蘇州醫(yī)工所碩士研究生.主要研究方向為腦電腦磁源成像.

E-mail：liutingcumt@163.com

戴亞康(通訊作者) 男,1982年出生于江蘇常州,現(xiàn)為中科院蘇州醫(yī)工所研究員、博士生導師.主要研究方向為醫(yī)學影像處理.在國內(nèi)外發(fā)表學術(shù)論文30余篇.

E-mail:daiyk@sibet.ac.cn

楊瑩雪 女,1986年出生于山東,現(xiàn)為首都醫(yī)科大學宣武醫(yī)院醫(yī)師.主要研究方向為腦功能疾病的臨床及電生理研究.

E-mail:yyx19861213@163.com

王玉平(通訊作者) 男,1961年出生于河北,現(xiàn)為首都醫(yī)科大學宣武醫(yī)院神經(jīng)內(nèi)科主任、北京市癲癇診療中心主任、腦功能疾病調(diào)控治療北京市重點實驗室主任、北京市腦重大疾病研究院癲癇病研究副所長.主要從事腦功能疾病的臨床、電生理和基礎(chǔ)研究,對癲癇、睡眠障礙、運動障礙病、認知障礙等臨床問題有較深入的研究.在國內(nèi)外發(fā)表學術(shù)論文280余篇.

E-mail:wangyuping01@sina.cn

An MEG Inverse Solver by Imposition of Temporal Smoothness Constraint

LIU Ting1,2,3,DAI Ya-kang1,YANG Ying-xue4,5,WANG Yu-ping4,5

(1.SuzhouInstituteofBiomedicalEngineeringandTechnology,ChineseAcademyofSciences,Suzhou,Jiangsu215163,China; 2.ChangchunInstituteofOptics,FineMechanicsandPhysics,ChineseAcademyofSciences,Changchun,Jilin130033,China; 3.UniversityofChineseAcademyofSciences,Beijing100049,China; 4.DepartmentofNeurology,XuanwuHospital,CapitalMedicalUniversity,Beijing100053,China; 5.BeijingKeyLaboratoryofNeuromodulation,Beijing100053,China)

The magnetoencephalography (MEG) inverse problem refers to the reconstruction of the neural activity of the brain from MEG measurements.A method to solve the MEG inverse problem employing temporal smoothness constraint is proposed under the assumption that time course of the source is smooth in time.Specifically,the temporal smoothness of the source was ensured by imposing a roughness penalty in the minimum norm estimate (MNE) data fitting criterion in the form of dual-parameter regularization.To select two tuning parameters,the generalized cross-validation criterion (GCV) was used.The inverse solutions were obtained as the linear combination of the one-parameter regularized solutions.We evaluated the proposed method by a synthetic example and a real data example.Compared with MNE,the proposed method can get smaller overall mean squared error (MSE) and smaller curvature variability.Moreover,the proposed method can reconstruct the shape of the time course of source better.

magnetoencephalography (MEG) time course;inverse problem;two-parameter regularization;temporal smoothness

2015-05-05;

2015-07-15;責任編輯:覃懷銀

中國科學院百人計劃基金項目;國家高技術(shù)研究發(fā)展計劃(863計劃)( No.2015AA020514);國家自然科學基金(No.61301042 );腦功能疾病調(diào)控治療北京市重點實驗室開放課題;江蘇省自然科學基金(No.BK2012189);蘇州市醫(yī)療器械與新醫(yī)藥專項基金(No.ZXY201426);中法“蔡元培”項目(No.201404490123)

TP301

A

0372-2112 (2016)12-2823-06

??學報URL:http://www.ejournal.org.cn

10.3969/j.issn.0372-2112.2016.12.002