基于腦磁信號的頭部三維能量分布模型研究*

李 東 柳 淵 孫 偉 白 玫*

基于腦磁信號的頭部三維能量分布模型研究*

李 東①柳 淵①孫 偉②白 玫①*

目的：通過分析處理腦磁數據，顯示全腦三維能量分布，以期準確的反映腦部的神經活動情況。方法：選擇小波函數對腦磁信號進行分析得到頻域數據，將能量信息映射到與之對應的傳感器位置，在傳感器之間進行適當的插值，創建三維坐標下的補片，顯示全腦三維能量分布。為驗證方法的有效性，采集手指運動試驗的腦磁數據，選取特定頻率段，分析討論不同時間下的地形圖。結果：能量分布圖反映的信息與運動事件相關同步研究的描述一致，即對側初級運動區(M1)和雙側輔助運動區(SMA)會出現高能量信號，神經活動活躍。結論：該方法可以全面準確地反映大腦的神經活動情況，為分析大腦神經活動、確定神經活動源發位置提供了有價值的參考方法。

腦磁圖；小波變換；地形圖

[First-author’s address] Department of Medical Engineering, Xuanwu Hospital, Capital Medical University, Beijing 100053, China.

人體大腦產生的磁場包括兩種，恒定磁場和交變磁場。恒定磁場由大腦內磁性物質產生，交變磁場則由神經細胞活動產生[1]。腦磁圖測量大腦神經細胞內電流產生的相關磁場，由腦內神經椎體細胞的電生理活動所產生。磁信號穿透腦組織和顱骨有很少的衰減或沒有衰減，與測量淺層細胞外電流靠容積傳導的腦電圖相比，腦磁圖具有較高的空間分辨率[2]。腦磁信號比腦電信號對于大腦皮層表面的活動更加敏感[3]。在腦神經活動磁源定位的空間分辨率上，腦磁圖優于腦電圖[4-5]。腦磁圖與腦電圖相比有著噪聲干擾小、信號采集優質及時間分辨率高等特點。近年來，腦磁數據被應用于神經外科術前定位和腦功能損傷評價等方面[6]。因此，對較為精準的腦磁信號進行分析，希望可以得到更加理想的處理結果。

1 腦磁信號處理方案

腦磁信號是非平穩信號，適合使用同時表達信號時域和頻域特性的時頻分析方法來處理。小波變換是典型的時頻分析方法，其用小波函數移位和縮放得到一系列小波表示的信號。在信號分解過程中，小波函數不斷平移、伸縮，使信號在各個頻率層面上得到分解。頻率從低到高，小波函數窗口會逐漸縮小。其分析結果在低頻段擁有高頻率分辨率和低時間分辨率，而在高頻段有低頻率分辨率和高時間分辨率。然而，小波函數在相當程度上影響變換的精準度和有效性。基于若干次先行試驗，選取Morlet小波，其在時頻分析處理腦電圖(electroencephalogram，EEG)信號時經常被使用[7]。

采取的處理方案步驟：①采集原始腦磁數據；②信號的預處理，為去除接觸噪聲、基線漂移和工頻干擾，對原始腦磁信號進行小波閾值法濾除噪聲；③選取合適的小波函數及窗函數參數，對腦磁各導數據進行小波分析；④數據插值，計算功率分布；⑤生成能量分布的三維偽彩色圖。而后還將引用1例手指運動試驗數據對該方法有效性進行驗證[8]。

2 三維頭模型的建立

2.1 腦磁圖系統傳感器位置排布

所用數據是由顱外102信道探測陣列所采集。每組傳感器測量三種數據，磁場強度、軸向梯度和曲面梯度。對于腦神經活動的淺表層圖像繪制，曲面梯度數據要優于軸向梯度數據。而對于定位大腦深部的場源，場強數據優于梯度數據[9-10]。故使用磁場強度信號(如圖1所示)。

圖1 顱外102組傳感器三維分布側面圖

2.2 數據插值方法

通過小波變換分析得到102組頻率域數據。以其顯示全腦能量分布略顯單薄，因此進行適當的插值。基于大數據量，將三維表示化簡為類似二維圖像插值的處理方法。傳感器在垂直軸方向分為7層排列，在相鄰兩層中各取相鄰兩點，4點上下相鄰，組成一個柵格。如圖2所示，A、B、C、D四點是相鄰的4組傳感器，其坐標分別為(xα，yα，zα)，α=A，B，C，D。X軸和Y軸方向插值數分別是n和m。

圖2 數據插值圖

位于4邊上的插值點可直接根據4點坐標計算得出。而后，線性插值計算出AB間n個插值點Sj，同理可得圖中的Tj，Pi，Qi，j=1，2，…n；i=1，2，…m。由此在柵格中形成一m*n的矩陣點集。計算出Pi與Qi之間第j個插值點Eij(公式1，公式2)。

式中j=1，2，…n。

式中i=1，2，…m，j=1，2，…n。

參數m和n的選擇需適當。過小則插值不會發揮其填充數據缺失的作用，若過大則會導致程序運行緩慢。考慮插值點數與運行代價的關系，最終選擇m=n=10[10]。插值后分布如圖3所示。

圖3 插值點俯視圖

在空間位置上越鄰近的點越有可能具有相似的特征值，這是空間插值技術中最基本的理論假設[11]。各種插值方法其共性是插值點幅值相當于對鄰近的原始數據加權求和。據此建立第m個插值點能量幅值計算模型(公式3)。

式中fi是原始數據點幅值；d表示插值點與原始數據點之間的距離；c(d)是與距離相關的權重函數；計算鄰近n個原始數據點的影響。

某點的能量強度對周邊幅值的影響呈多種形態分布，這里取高斯分布，原始數據的能量幅值A可作(公式4)：

式中x為插值點與數據點間的距離。則插值點能量和相距d的數據點間的關系為(公式5)：

然而，某一插值點的能量值是原始數據點對其影響的總和。遠距離的數據點對插值點的能量貢獻甚微，可忽略，故取最臨近的4個原始數據點計算。

使用MATLAB編寫程序。創建補片，一個補片對象是由其頂點確定的多邊形。如果坐標數據是不能定義封閉的多邊形，將自動使多邊形封閉。然而，如果單個補片面的邊緣相互交叉，得到的面可能不會完全填充。此時，應將面分解為更小的多邊形；補片色譜自動取偽彩色。

3 臨床實驗數據的處理

3.1 手指運動實驗

手指活動激活對側中央區和頂區皮層，皮層活動的10～25 Hz頻率成分首先在-0.5～0.3 s時間段表現為功率抑制，接著有強烈的功率增強活動顯示，最顯著的功率增強集中在0.4～0.7 s區間[12]。事件相關功能磁共振成像顯示初級運動區(M1)為運動執行區，手指運動對側M1區明顯激活，同側則無明顯激活區；輔助運動區(supplementary motor areas，SMA)、運動前區(premotor areas，PMA)、基底節及小腦皮質均為雙側激活[13]。基于上述研究，對三維能量分布的頭模型進行實驗驗證。本研究實驗數據來自首都醫科大學宣武醫院腦磁圖室，為手指運動實驗數據。被試手指為右利手，無手部骨關節疾病，無金屬義肢、義齒，精神狀態良好。實驗中要求受試者保持安靜，閉眼，盡量避免眼動、咬牙及皺眉等動作；左、右手食指分別平放在兩個響應按鍵上；給予單耳純音刺激，隨機給予左、右耳刺激；根據聽到的刺激，同側食指做出按壓鍵反應，其余手指不動。所用為右手運動數據。

3.2 結果分析

腦磁信號呈現特征節律波，可在枕部和感覺運動區記錄到α節律和β節律，其節律表現出事件相關變化，若對腦磁信號進行頻率域分析，可進一步了解大腦特定的活動機制[14]。試驗數據采樣頻率1000 Hz，對刺激前0.1 s到發出指令后2.4 s給予關注。

α節律取10 Hz為代表頻率，時間取第0.05 s、0.1 s、0.2 s、0.3 s、0.4 s、0.5 s、0.6 s、0.7 s、0.8 s、0.9 s和1.0 s。即指令發出前0.05 s到發出后0.9 s的數據。兩圖呈現能量的變化如圖4、圖5所示。

在0.05 s高能量信號分布在左側中央溝、中央前回和中央后回區域，可理解是M1對應的部分。第0.1～0.6 s扣帶回前部出現高能量分布，能量最高值出現在0.4～0.5 s，可理解其對應為雙側SMA。第0.7～1.0 s雙側SMA能量降低。指令發出之后，即0.1 s開始，能量開始逐漸增加，在0.4～0.5 s之間有顯著的增加，0.5 s開始能量下降，持續至第0.6 s。第0.6 s開始顯著下降。到第0.8 s，能量狀態恢復至0.1 s之前。整個過程10 Hz節律活動最為強烈的時間是在指令發出后的第0.4 s。圖5 f給出了10 Hz能量最值隨時間變化的情況。大致可將其分為三個階段：0～0.4 s，能量稍有升高；0.4～0.6 s，能量顯著升高；0.6 s以后逐漸恢復[15]。

圖4 10 Hz信號0.05～0.5 s時間點能量三維分布圖

圖5 10 Hz信號0.6～1.0 s時間點能量三維分布圖

β節律取25 Hz，時間取第0.05 s、0.1 s、0.2 s、0.3 s、0.4 s和0.5 s。即指令發出前0.1 s到發出后0.4 s的數據。圖6呈現β節律段25 Hz成分信號能量的全腦三維分布隨時間變化的情況。同樣開始時，高能量集中在對側M1。之后的0.3 s內，高能量信號出現在雙側SMA，且能量值逐漸增大。到指令發出后的0.4 s，雙側SMA高能量信號減弱消失。0.1 s開始能量逐漸升高，這一趨勢直至第0.5 s結束。

圖6 25 Hz信號0.05～0.5s時間點能量三維分布圖

文獻研究顯示，對于利手的運動事件相關變化，主要體現在對側M1，雙側SMA及小腦的部分[15]。這與前述的分析討論比較一致。

4 結論

皮層神經活動在時頻域的能量分布情況可體現事件相關的皮層神經活動特征。本實驗結果與有關文獻對于手指運動的事件相關同步的若干研究情況相吻合[13-15]。采用三維考察腦部神經活動的方法，可觀察三維空間中腦磁信號不同時間、不同頻率下的能量分布地形圖，以分析大腦的神經活動情況[16]。為分析大腦神經活動源發位置提供了較為有價值的觀察方法。由于條件所限，本研究仍然存在某些問題，如能量分布圖有明顯的插值痕跡，插值方法有待改進。今后還需細化插值方式，提高程序運行速度，以設計出可視化的界面，提供更加細致的顯示形式。

[1]韓建坤，康文巧.腦磁圖簡介[J].醫療裝備，2004，17(8)：15-17.

[2]Ballief S，Mosher J C，Leahy R M.Electromagnetic brain mapping[J].IEEE Signal Processing Magazine，2001，18(6)：14-30.

[3]王建軍.腦圖和腦電圖在癲癇中的應用[J].臨床神經電生理學雜志，2005，14(3)：180-183.

[4]Cuffin BN.Effects of local variations in skull and scalp thickness on EEG's and MEG's[J].IEEE Transactions on Biomedical Engineering，1993，40(1)：42-48.

[5]Nicol'as VE，Carlos HM，Michael W，et al.Effect of head shape variations among individuals on the EEG/MEG forward and inverse problems[J].IEEE Transactions on Biomedical Engineering，2009，56(3)：578-596.

[6]吳婷，劉宏毅.腦磁圖在神經外科術前定位的研究進展[J].現代電生理學雜志，2011，18(2)：102-105.

[7]張萍，黃慶.腦磁圖表現對腦卒中患者腦功能損傷程度和區域的評估價值[J].中國臨床康復，2004，8(34)：7650-7651.

[8]吳國材，李梅，吳南，等.腦磁圖、視頻腦電圖在癲癇術前評估中的應用[J].第三軍醫大學學報，2009，31(11)：1080-1083.

[9]崔驪，李向東，云慶輝.心腦電圖機測量結果的不確定度分析評定[J].中國醫學裝備，2013，10(6)：26-27.

[10]柳淵，孫偉，嚴漢民，等.癲癇患者腦磁圖的高頻振蕩量化算法研究[J].中國醫學裝備，2012，9(11)：1-3.

[11]Zygierewicz J，Sieluzycki C，Konig R，et al.Eventrelated desynchronization and synchronization in MEG：Framework for analysis and illustrative datasets related to discrimination of frequency-modulated tones[J].J Neurosci Methods，2008，168(1)：239-247.

[12]鄔倫，劉瑜，張晶.地理信息系統原理方法和應用[M].北京：科學出版社，2005：178-192.

[13]侯文生，李衛娜.Jiang Yingtao，等.手指活動相關腦磁圖的時-頻域特征分析[J].儀器儀表學報，2009，30(10)：2093-2098.

[14]李恩中，韓瓔，翁旭初，等.事件相關功能MRI在隨意運動腦功能活動區協同作用研究中的應用[J].中國放射學雜志，2002，36(5)：397-401.

[15]Pollok B，Gross J，Schnitzler A.How the brain controls repetitive finger movements[J].Journal of Physiology-Paris，2006，99(1)：8-13.

[16]于薇，林沖宇，臧玉峰，等.利手和非利手隨意運動的全腦功能磁共振成像[J].中國放射學雜志，2003，37(5)：402-405.

Research on head three-dimensional energy distribution model based on brain magnetoencephalography signal processing

/
LI Dong, LIU Yuan, SUN-Wei, et al// China Medical Equipment,2014,11(10):8-11.

Objective: A three-dimensional energy distribution over the whole brain will be displayed by analyzing the brain magnetic data, so as to reflect the neurological activity of the brain accurately. Methods: The brain magnetic signal was processed with wavelet function to get frequency domain data. The energy information is mapped to the corresponding sensors’position. In order to show the energy distribution, the author interpolated between sensors, created three-dimensional patches. For the purpose of validating, MEG data from a finger exercise testing subject was analyzed. Moreover, data of particular frequency was discussed also. Results: What it shows is consistent with description on event related synchronization. Actually the movements activated contralateral primary motor cortex (MI), bilateral supplementary motor area (SMA). Conclusion: Because of reflecting neurological activity accurately, it is said this method provides a valuable reference to analyze brain activity and locate neural activity.

Magnetoencephalography; Wavelet transform; Topographic map

1672-8270(2014)10-008-04

R831

A

10.3969/J.ISSN.1672-8270.2014.10.003

2014-06-26

國家自然科學基金(81171229)“應用腦電磁信號頻段分解對癲癇病人腦皮層中央區振蕩規律的研究”

①首都醫科大學宣武醫院醫學工程處 北京 100053

②首都醫科大學宣武醫院神經內科 北京 100053

*通訊作者：jswei65@163.com

李東，女，(1987- )，碩士，助理工程師。首都醫科大學宣武醫院醫學工程處，研究方向：醫學信息分析與處理。