李明彩，張業(yè)宏*，秦 鑫

（1. 新鄉(xiāng)醫(yī)學院醫(yī)學工程學院，新鄉(xiāng) 453003；2. 河南省神經(jīng)傳感與控制工程技術(shù)研究中心，新鄉(xiāng) 453003）

0 引言

人類腦計劃最早由美國于2013 年啟動，近年來腦科學的巨大潛力持續(xù)受到各國政府、科學界、商界的高度重視并競相啟動腦計劃任務(wù)。我國更是把腦科學與類腦研究列為2030 重大項目及十四五規(guī)劃國家重大科技項目。腦科學面對的是人類社會的基礎(chǔ)科學問題之一，也是人類亟待深入探索的領(lǐng)域。大腦是由上百億神經(jīng)元和多級神經(jīng)系統(tǒng)組成的復雜信息網(wǎng)絡(luò)，網(wǎng)絡(luò)中靠腦電進行信息的傳播，有效采集腦電并探討其成分與特定生理活動相關(guān)性是研究大腦信息處理機制的重要環(huán)節(jié)。

植入式腦電采集技術(shù)能將微電極深入生物體大腦皮層，記錄到細胞外神經(jīng)元水平的電信號，局部場電位（local field potential, LFP）就是電極尖端區(qū)域附近興奮性和抑制性神經(jīng)元群體電位信號總和，反映了局部網(wǎng)絡(luò)中各神經(jīng)元的協(xié)同作用。因其屬于有損采集，實驗中主要以嚙齒類動物如小鼠為主，參考人類腦電低頻信號的頻段分類，常將LFP 低頻部分劃分為delta波段（0.5～4 Hz）、theta 波段（4～8 Hz）、alpha波段（8～13 Hz）、beta 波段（13～30 Hz）。已有研究表明，delta 腦電波段的震蕩與小鼠清醒狀態(tài)下晶須桶狀皮層局部場電位的活動被呼吸鎖相，delta 波與呼吸活動存在相關(guān)性［1］。陸續(xù)有學者采用獨立成分分析、小波變換、經(jīng)驗模態(tài)分解等腦電波分離技術(shù)分離出delta波段信號，與呼吸信號進行相關(guān)系數(shù)、互功率譜等計算［2-4］。曾麗萍［5］以腦電波分離和濾波技術(shù)為基礎(chǔ)計算小鼠平均呼吸發(fā)生率和δ波段的積差相關(guān)系數(shù)。由于腦電成分的復雜性、頻率混疊等問題，局部場電位信號具有典型的非線性特性，通過簡單的計算相關(guān)系數(shù)或互功率譜，總體相關(guān)性并不顯著。吳孫勇等［6］采用獨立成分分析對小鼠腦電波進行分離得到delta波，進而進行希爾伯特-黃變換得到信號的瞬時頻率信息，將瞬時頻率超出閾值的異常時刻與呼吸信號波峰時刻建立相關(guān)性分析，得出呼吸信號與delta波段顯著相關(guān)的結(jié)論。本文則在不進行弱腦電波信號分離的前提下，采用小波變換對腦電和呼吸信號進行時頻域分解，進而通過時頻互信息法來定量分析腦電波特定波段與呼吸過程的相關(guān)性。

1 實驗數(shù)據(jù)描述和分析

實驗數(shù)據(jù)來源于2014 年全國研究生數(shù)學建模競賽，小鼠因體型小植入式腦電采集到通道數(shù)相對較少，五個植入式電極以0.25 mm為間隔、以1 kHz為采樣頻率并排記錄得到小鼠在睡眠狀態(tài) 下1～5 通 道（LFPchl1、LFPchl2、LFPchl3、LFPchl4、LFPchl5）視覺皮層局部電位，由于嚙齒類小鼠僅需考慮鼻呼吸，利用鼻子附近外置熱敏電極間接記錄第6通道，得到呼吸曲線，其幅值沒有生理意義，而局部極大波峰值和極小波谷值時間點標記了呼氣周期和吸氣周期的開始，利用MatlabR2021a軟件繪制了小鼠睡眠狀態(tài)下原始呼吸曲線及局部極值，如圖1所示，實驗記錄20 s時長呼吸信號，通過統(tǒng)計局部極值計算峰峰時間間隔可知平均周期為0.5695 s，平均頻率為1.7643 Hz，用快速傅里葉變換計算呼吸信號功率譜，分析結(jié)果如圖2所示，其第一個低頻峰值約1.8 Hz，對應(yīng)于呼吸過程節(jié)律變化最具周期性的部分，時域與頻域周期分析具有一致性。

圖1 睡眠狀態(tài)下小鼠呼吸曲線及局部極值

圖2 原始呼吸信號功率譜

為分析視覺皮層局部場電位腦電信號的周期節(jié)律，分別采用自相關(guān)函數(shù)和功率譜分析的方法從時域頻域進行討論，自相關(guān)函數(shù)在一定程度上刻畫了每個序列在不同相位差的情況下與自身相似性，相位差為0 時達到最大，圖3 展示了LPF 中1～4 通道局部場電位信號自相關(guān)函數(shù)運算結(jié)果，可以觀察到睡眠狀態(tài)下小鼠每個通道的LFP 信號都具有相似的起伏規(guī)律，體現(xiàn)較強的自相關(guān)性，表明其隱含有周期節(jié)律。圖4展示了通道1 信號的功率譜，如圖4 可知，原始LFP 信號頻率主要分布在0～4 Hz 之間，主要成分為delta 波，但從中不能定位出與呼吸峰值區(qū)域?qū)?yīng)的成分。由于LFP 信號低頻低幅值、非平穩(wěn)、隨機性強，關(guān)注的信號會被較強的背景信號掩蓋，只有分離出與呼吸活動相關(guān)聯(lián)的腦電信號后再計算信號的時域相關(guān)系數(shù)才是有意義的［7］。相較于時域特征，LFP 信號頻域特征更突出，考慮通過小波變換將信號轉(zhuǎn)換到時頻域并采用互信息法探討二者的相關(guān)性，著重關(guān)注不同通道信號中相應(yīng)頻率成分的信號是否具有相似的波動特征。

圖3 睡眠狀態(tài)下LFP1～4通道腦電波自相關(guān)函數(shù)圖

圖4 初級視覺皮層局部場電位功率譜（LFPchl1）

2 呼吸與腦電信號相關(guān)性分析

2.1 小波變換

為研究不同通道信號特定頻率成分的相關(guān)性，通過前文分析，采用小波變換將相應(yīng)腦電和呼吸信號轉(zhuǎn)換到時頻域以計算特定頻率的強度特征，小波變換通過引入平移因子τ和伸縮不變性尺度參數(shù)a實現(xiàn)可變窗口，如式（1）所示，靈活選取恰當?shù)男〔ɑ瘮?shù)，用一簇函數(shù)去逼近原信號，其多分辨率特性可實現(xiàn)時頻窗口隨信號特征的自適應(yīng)改變，小波變換使其具有局部分析和細化功能［8］。因無需重構(gòu)變換而采用具有對稱性的Morlet小波基函數(shù)，如式（2）所示。

2.2 互信息

互信息（mutual information, MI）的概念來自于概率論與信息論，反映了兩個隨機變量間線性或非線性關(guān)聯(lián)的程度?；バ畔⒌娜≈凳欠秦摰模瑸?表示兩個隨機變量相互獨；若歸一化后達到1，則表明二者關(guān)聯(lián)度最大［9］。假定Xi，Xj均為離散隨機變量，Bins表示隨機變量分布范圍內(nèi)劃分的區(qū)域數(shù)量（b=1,2,…Bins），采用直方圖法統(tǒng)計隨機變量在各區(qū)域的分布，可估計出其概率密度函數(shù)分別為p（Xi）和p（Xj），進而可定義隨機變量的信息熵H和二者互信息計算，分別如式（3）和式（4）所示，p（Xi,b,Xj,b）為其聯(lián)合概率密度函數(shù)，表示兩個事件同時發(fā)生的概率，表現(xiàn)為Bins*Bins大小的數(shù)據(jù)矩陣。

2.3 時頻互信息法計算相關(guān)性算法步驟

第一步：以圖1中呼吸信號極大值時刻為起點，分別對1～5 通道腦電和第6 通道呼吸信號進行分段處理，每個通道的信號都被分成S段長度為N點的時間序列子段，每個子段信號仍按6個通道組合在一起并歸一化（歸一化的原因是腦電信號和呼吸信號幅值差異較大）為6*N大小的矩陣X=(x1,x2,x3,x4,x5,x6)，xi（i=1, 2, …, 6）均為行向量，x1～x5為腦電信號，x6為呼吸信號；

第二步：根據(jù)式（1）、式（2）對xi進行指定頻帶［flow，fhigh］時頻域小波分解，頻率分辨步長為fstep，得到每段信號的時頻能量譜圖，大小為（fhigh-flow）/fstep*N，沿頻率軸方向得到信號平均能量波動向量Fi(i= 1,2,…,6)；

第三步：根據(jù)式（3）、式（4）分別求取每個離散隨機變量Fi的概率密度函數(shù)p（Fi）、聯(lián)合概率密度函數(shù)p(Fi,b,Fj,b)及互信息I(Fi,Fj)，其中i=1,2,…,5，j=6，b=1,2,…,Bins；

第四步：第一步S段信號組合重復執(zhí)行第二步和第三步，對所有分段得到的互信息值I(Fi,Fj)取均值并歸一化Imeans最終求得5通道腦電與呼吸信號互信息值。

2.4 算法結(jié)果及分析

根據(jù)2.3 節(jié)步驟進行算法實現(xiàn)，經(jīng)過對分段信號長度N、頻率分辨步長fstep、概率密度函數(shù)分布區(qū)域數(shù)量Bins等參數(shù)進行調(diào)整并反復測試。由于每個呼吸周期時長均有差異，N可在500～600 之間取值。根據(jù)前文對信號的功率譜分析fstep取0.1 Hz 在頻率分辨率和時間復雜度上較為合適，Bins按照概率密度函數(shù)常用要求取40～60計算。圖5為第二步的計算結(jié)果，以通道2第二段腦電信號時頻能量譜示例，與前文腦電頻域分析結(jié)果一致；圖6為第三步的計算結(jié)果，以通道2 第二段s=2 腦電信號與第二段呼吸信號聯(lián)合概率密度分布示例，其聯(lián)合分布具有顯著規(guī)律。表1展示了第四步時頻互信息相關(guān)性分析的計算結(jié)果。從表1可知，當把腦電信號時頻分析的頻帶設(shè)在波段［0.5，4］Hz、呼吸信號時頻分析的頻帶設(shè)在波段［0.5，2］Hz 時，每個通道腦電信號與呼吸信號的互信息值均明顯高于其他波段，二者在相應(yīng)頻帶的成分具有最相似的時頻波動特征，其平均互信息值歸一化后達到0.7672，因此得出結(jié)論：呼吸節(jié)律與0.5～4 Hz 范圍的delta波顯著相關(guān)。

表1 基于時頻互信息法的腦電與呼吸信號相關(guān)性分析結(jié)果（N=510，fstep=0.1 Hz，Bins=60）

圖5 通道2腦電信號s=2時頻能量譜圖

圖6 通道2腦電s=2和呼吸信號s=2聯(lián)合概率密度分布

3 結(jié)語

為討論小鼠睡眠狀態(tài)下植入式電極采集到的初級視覺皮層局部場電位腦電信號和呼吸生理活動的相關(guān)性，將分段后的兩類信號進行小波變換以求取信號在時頻域中的平均能量波動，進而以時頻域平均能量波動為離散隨機變量計算其概率分布密度和聯(lián)合概率分布密度，繼而求得兩種分布特征之間的互信息值，最終得到delta 波和呼吸信號時頻互信息平均值為0.7672，因此認為呼吸信號和delta 波存在顯著性相關(guān)。本文為探索腦電信號成分與特定生理活動相關(guān)性提供了有價值的思路和方法。