快速有效地提取和比較不同神經信號中所含信息的相同和區別之處是研究者關注的問題。本文介紹了一種基于小波變換的時頻建模方法。該方法用少量的半橢球模型來表征神經信號在時域和頻域上的變化,克服了傳統時頻分析中背景干擾大、參數多的缺陷。將該方法應用于嗅球場電位的研究,與支持向量機(SVM)算法相結合,可初步實現氣味的分類。
引用本文： 董琪, 胡靚, 莊柳靜, 周俊, 王平. 基于小波變換的時頻圖建模及其在嗅覺場電位分析中的應用. 生物醫學工程學雜志, 2014, 31(3): 481-486. doi: 10.7507/1001-5515.20140089 復制
引言
近年來，隨著腦-機接口技術的發展，低頻神經信號(＜200 Hz)，包括非侵入式檢測(如腦電圖，electroencephalogram，EEG)[1] 和侵入式檢測(如局部場電位，local field potential，LFP)[2]，得到越來越多研究者的重視。和單神經元放電(spike)信號相比，這些低頻神經信號反映了一定區域內所有神經元電生理活動的總體特征，因此具有更高的穩定性，并被廣泛應用于神經解碼、信號分析、疾病檢測等領域 [3-8]。
時頻聯合分析是上世紀90年代發展起來的神經信號處理方法。該方法能夠表征神經信號的頻率成分組成以及相關頻率成分的能量分布隨時間的變化，從而克服了傳統時域和頻域分析不能完整準確地表征信號的缺陷[1]，在各種神經信號研究中得到了大量的應用[5-7, 9]。在EEG和LFP的分析中，主要的難點在于：(1)與事件相關的信息往往被淹沒在強大的背景噪聲中，很難從原始信號中有效提取有用特征[3]；(2)基于短時傅里葉變換等手段的傳統時頻分析方法會產生大量的參數，不利于統計歸納。
本文介紹了一種新的低頻神經信號的處理方法。該方法將小波時頻圖和橢球方程建模方法有機地結合起來，用時頻圖表征神經信號的時頻特征，并通過建模的方式從雜亂的背景噪聲中有效提取其中的有用信息，大大減少了時頻圖原有的參數量，在一定程度上克服了傳統時頻分析中存在的缺陷。將該方法應用于大鼠嗅球LFP的研究，可提取與氣味刺激相關的事件特征，與支持向量機(support vector machine,SVM)算法相結合，初步實現了三種氣味刺激的分類。
1 分析方法
該方法的核心思想是用盡可能少的方程參數來表征原始時頻圖，以達到去除噪聲干擾，提取有用信息供進一步分析的目的。主要包括兩個步驟：① 構建神經信號基于小波變換的功率譜時頻圖；② 使用半橢球方程對時頻圖進行建模。
1.1 基于小波變換的功率譜時頻圖
小波變換是以某些特殊函數為基礎，將數據過程或數據系列變換為級數系列，以發現它的類似頻譜的特征，從而實現數據處理。這種方法是由法國工程師Morlet 在1974 年首先提出，經過之后數十年的研究，該方法建立了完整的數學化體系和扎實的理論基礎，并被廣泛應用于各種神經生理信號的分析[1-2, 10-13]。函數f(x)的連續小波變換可定義為
| ${{W}_{f}}\left( a,b \right)=\frac{1}{\sqrt{\left| a \right|}}\int_{-\infty }^{+\infty }{f\left( x \right)\bar{\Psi }}\left( \frac{x-b}{a} \right)dx,$ |
式中Ψ(x)是母小波，(x)是Ψ(x)的共軛形式，參數a 和b 分別被稱為尺度因子和平移因子。由式(1)可知，連續小波變換是空間(時間)和尺度的局部變換，因而能有效地從信號中提取信息。通過伸縮和平移等運算功能可對函數或信號進行多尺度的細化分析，解決了傅里葉變換不能解決的許多困難。
由連續小波變換得到的小波功率譜可定義為
| $TF\left( a,b \right)={{\left| {{W}_{f}}\left( a,b \right) \right|}^{2}}$ |
式中TF(a,b)是函數f(x)在尺度a和位移b時所具有的能量。嚴格地說，連續小波變換提供了時間尺度分析，而不是時間頻率分析。然而，使用適當的尺度頻率變換，依然可以得到一個近似于時間頻率分析的工具。與短時傅里葉變換相比，小波時頻圖的窗函數具有可變的特性。對于較高的頻率，時間窗會變小，使時間分析精度提高，而對于較低的頻率，時間窗會變大，保留了完整的低頻信息。
小波變換的關鍵是要選擇合適的母小波[1]。復Morlet小波是一種使用廣泛的基本小波,已被眾多研究證明，適用于電生理信號的時頻分析[1-3, 10-11, 13]。復Morlet小波w(t,f)對于某頻率的表達式為
| $w\left( t,f \right)=Aexp(2i\pi ft)exp({{t}^{2}}/2{{\sigma }_{i}}^{2})$ |
式中A為歸一化因子，其目的是使小波能量為1；f與σt(時域分辨率)是需要適當選擇的參數，兩者的乘積決定了小波的周期數，同一復Morlet小波家族具有恒定的f與σt乘積，本研究所使用的小波家族定義為2πfσt=7(見圖 1a)。復Morlet小波w(t,f)是一種復值調制的高斯(Gaussian)函數，在時域和頻域上都具有高斯分布，其傅里葉頻譜的標準差滿足
| ${{\sigma }_{f}}=\frac{1}{2\pi {{\sigma }_{t}}},$ |
特別的，在本研究中，
| ${{\sigma }_{f}}=\frac{1}{2\pi {{\sigma }_{t}}}=\frac{f}{7},$ |
因此，不同頻率f所對應的σt和σf不相同，即在整個時頻平面上具有可變的時頻分辨率：在高頻區能提供高的時間分辨率，在低頻區能提供高的頻率分辨率(見圖 1b)。
圖1
復Morlet小波的特性
(a)復Morlet小波的實部(公式3)，左：f/
(a) real part of equation 3,with
由于神經信號在低頻段具有的能量比高頻段大很多，為了在同一比例下觀察場電位信號在低頻段和高頻段能量隨時間的變化，由連續小波變換獲得的信號功率時頻圖需要標準化
| ${{z}_{ft}}=\frac{{{c}_{ft}}-{{m}_{f}}}{{{\delta }_{f}}},$ |
式中cft為f(t)在頻率f及時間t時的小波功率譜，mf和δf分別為刺激前1 s內信號在頻率f的小波功率譜的平均值和標準差。
圖 2為一段人造信號的原始波形及其基于復Morlet小波變換的時頻圖。該人造信號(見圖 2a)模擬了神經信號在各頻率段的成分組成，具有1個5 Hz的基礎振蕩，并在0.5~0.7 s和1.5 ~1.7 s范圍內分別具有一個β段(15~40 Hz)和σ段(40~90 Hz)的突變。圖 2(b)為該信號的原始時頻圖，由于低頻段基礎振蕩的能量大，使得高頻段的突變不明顯，經公式(6)歸一化之后，能量的突變點被凸顯出來，而恒定的能量波段被弱化。
圖2
人造信號及其小波時頻圖
(a)人造信號原始波形；(b)人造信號原始小波時頻圖；(c)經過歸一化之后的小波時頻圖
Figure2. An artificial signal and its wavelet time-frequency maps(a) the waveform of an artificial signal； (b) the raw wavelet time-frequency map； (c) the normalized wavelet time-frequency map
1.2 基于半橢球方程的時頻圖建模
原始時頻圖包含的數據量非常大，很難直接提取信號的時頻特征。在本研究中，采用了建模的方法進行信息提取和降緯處理。該方法的核心思想是利用事先預定的參數方程φ對原始時頻圖進行擬合，主要包含時頻圖窗口化、求和與擬合三個步驟。
首先將原始時頻圖分解成一系列重疊的窗口，以確定需要建模的區域。在經行窗口化時，窗口的尺寸必須滿足三點要求：① 由于原始時頻圖是通過連續小波變換獲得，在進行窗口化時必須使窗口的尺寸能夠自動適應小波變換的時頻分辨率，即L/σt=H/σf。其中L和H分別為窗口的長度和高度，σt和σf分別為小波變換的時域標準差和頻域標準差。 ② 要選擇合適的窗口寬度，在充分體現頻率特征的基礎上凸顯能量的突變，一般選擇的窗口需要包含當前頻率3~5 個振蕩周期。計算中，我們選取的值為4。 ③ 為了使位于不同位置的窗口具有相同的評判標準，窗口的面積必須相同。
因此，對于特定的頻率f，窗口的長度L=4/f，即一個窗口包含當前頻率的4個振蕩周期，則L/σt=8π/7，由此換算出窗口的高度H=8πf/49。用這個值去計算原始時頻圖上每個點，就能獲得一系列面積為3249π的重疊窗口。然后逐一計算每個窗口中的能量總和，窗口總能量值體現了該窗口內場電位的振蕩強度。對窗口wi，總能量值Si計算公式如式7所示，然后找出所有大于預設閾值的窗口。
| ${{S}_{i}}=\sum\limits_{t,f\in {{w}_{i}}}{{{z}_{ft}}}$ |
最后，在選中的窗口中使用參數方程進行擬合。本研究選擇半橢球方程φ作為擬合方程并定義為
| $\varphi \left( f,t \right)=\left\{ \begin{align} & a\sqrt{1-v},0 ＜v＜1 \\ & 0,v>1 \\ \end{align} \right.$ |
其中v=(f/μf)2/lf2+(t/μt)/lt2，μf和μt為該橢球在時頻圖上的位置，a為橢球的極半徑(幅度)，lf和lt為半橢球的赤道半徑。圖 3為一個半橢球方程的示意圖。該方程共有5個參數，并且需要滿足以下條件：
(1)μf ＞ 0，μt ＞ 0，且使半橢球位于所選窗口的中心；
(2)0 ＜ lt ＜ L,0 ＜ lf ＜ H。
圖3
由5個參數進行表征的半橢球方程示意圖
Figure3.
Half-ellipsoid model with its five parameters
在擬合的過程中，要使建模誤差函數C達到最小，函數C定義為
| $C=\frac{1}{2}\sum\limits_{t,f\in {{w}_{i}}}{({{z}_{ft}}-\varphi \left( f,t \right)),}$ |
圖 4為采自麻醉大鼠在香芹酮刺激下嗅球的局部場電位β段信號的小波時頻圖以及建模后的結果。從圖中我們可以發現，用5個半橢球模型就能將原始時頻圖的主要能量變化從大量參數和背景噪聲中提取出來。
圖4
嗅球場電位信號β段時頻圖和半橢球方程建模結果
左：歸一化后的小波時頻圖；中：半橢球建模的2D表示；右：半橢球建模的3D表示
Figure4. Time-frequency map (left) of β range of signal from rat’s olfactory bulb and its modeling results in two (middle) and three (right) dimensions2 大鼠嗅球場電位記錄和氣味識別分析
氣味分子的識別是由嗅黏膜上的一千多種氣味受體完成的，一種氣味受體只能與特定分子結構的氣味相結合[15-17]。表達同種氣味受體的感覺神經元投射到嗅球表面的一個或幾個嗅小球中，因此特定的氣味能夠使嗅球產生特定時空特性的活動[18-22]。研究者們發現，氣味刺激前后嗅球僧帽細胞層神經電活動十分強烈，LFP振蕩會明顯增加，其幅值和波形的變化受到氣味種類的調控 [9, 23-25]。
2.1 實驗設計和信號記錄
實驗對象SD大鼠購于浙江省醫學科學院，質量是200~250 g。大鼠進行電極植入手術流程是在腹腔注射10% 的水合氯醛(4 mL/kg)麻醉后，移除單側嗅球上覆蓋的頭骨，利用油壓微推進器(Narishige,Japan)將自制的16通道直徑為65 μm鎳鎘合金微絲(AM system,USA)電極陣列，植入到嗅球僧帽細胞層，平均深度約為200~300 μm。由于嗅球僧帽細胞放電幅值較大，能提示電極已植入到目標細胞層，故我們選取氣味刺激時能記錄到細胞鋒電位信號的通道，分析該通道的低頻場電位信號。
使用濃度為10 mmol/L的異丁醇、苯甲醚和香芹酮三種氣味作為刺激物，單次刺激給予體積為13 mL、時間為5 s的氣味，刺激間隔為60 s。多通道信號采集裝置為商用OmniPlex數據采集系統(Plexon,USA)，具有20倍的前置放大器，16位的A/D轉換器，能以最高頻率為40 kHz，5 000倍模擬放大，實現256個通道實時數據采集，并帶有用于低頻場電位和高頻鋒電位的頻帶隔離可編程硬件濾波器。信號采集裝置通過網線將采集到的數據發送到計算機，利用Matlab編譯環境實現數據的處理和分析。實驗的整體示意圖如圖 5所示。
圖5
實驗整體示意圖
Figure5.
Schematic of the experiment
2.2 特征提取和氣味分類
考慮到了小波變換和半橢球方程建模時可能造成的邊界效應，我們將記錄的嗅球場電位信號分成三個獨立的頻段分別進行變換。頻率的步進分別為0.25 Hz(θ段、β段)和1 Hz(γ段)。1.5~18.5 Hz：用于研究θ段的信號(5~12 Hz)；9~37.5 Hz：用于研究β段的信號(12~30 Hz)；25~85 Hz：用于研究γ段的信號(30~75 Hz)。經過半橢球方程建模之后，原始信號可以用一系列的半橢球的參數來表征。我們定義了以下參數作為分類的特征值：F1，橢球的個數；F2，橢球的最大幅值; F3，橢球的幅值總和。
因此，每個頻段能獲得3個特征值，每段信號能獲得9個特征值。這些特征值經過歸一化之后，用于后續的基于SVM算法的氣味分類。 本文采用了SVM算法[26]對氣味信號行進分類識別。三種氣味共218組時頻圖向量(苯甲醚73組，異丁醇72組，香芹酮73組)，隨機選取一半作為樣本A，另一半作為樣本B，進行分類測試，最終識別結果如表 1所示。二重交叉檢驗得到的總體準確率為79.4%，其中苯甲醚、異丁醇和香芹酮的檢測準確率分別為82.2%、80.6%和75.3%，其準確度和特性性均高于在以往研究中基于短時傅里葉變換的分析和分類方法[9]。
表1
氣味分類結果
Table1.
Results of odor discrimination
3 結論
本文介紹了一種新型的低頻神經信號的處理方法。該方法將基于小波變換的時頻圖分析方法與半橢球建模方法有機結合起來，能夠快速有效地提取神經信號能量變化的有用信息，也可用于比較大量數據在相同與不同刺激下的相同與不同之處，克服了傳統時頻分析方法背景噪聲大、參數多的缺陷。該方法目前存在的主要問題是建模算法比較復雜，當待處理數據段增長時，對計算機硬件配置要求較高，數據處理時間大幅增加，有待進一步的優化。
將該方法用于麻醉大鼠嗅球場電位的分析中，可用于提取不同氣味刺激下，嗅球的電生理響應特征，通過模式識別的方法，初步實現不同氣味刺激的分類，其準確率和特異性均高于基于短時傅里葉變換的分析和分類方法。
引言
近年來，隨著腦-機接口技術的發展，低頻神經信號(＜200 Hz)，包括非侵入式檢測(如腦電圖，electroencephalogram，EEG)[1] 和侵入式檢測(如局部場電位，local field potential，LFP)[2]，得到越來越多研究者的重視。和單神經元放電(spike)信號相比，這些低頻神經信號反映了一定區域內所有神經元電生理活動的總體特征，因此具有更高的穩定性，并被廣泛應用于神經解碼、信號分析、疾病檢測等領域 [3-8]。
時頻聯合分析是上世紀90年代發展起來的神經信號處理方法。該方法能夠表征神經信號的頻率成分組成以及相關頻率成分的能量分布隨時間的變化，從而克服了傳統時域和頻域分析不能完整準確地表征信號的缺陷[1]，在各種神經信號研究中得到了大量的應用[5-7, 9]。在EEG和LFP的分析中，主要的難點在于：(1)與事件相關的信息往往被淹沒在強大的背景噪聲中，很難從原始信號中有效提取有用特征[3]；(2)基于短時傅里葉變換等手段的傳統時頻分析方法會產生大量的參數，不利于統計歸納。
本文介紹了一種新的低頻神經信號的處理方法。該方法將小波時頻圖和橢球方程建模方法有機地結合起來，用時頻圖表征神經信號的時頻特征，并通過建模的方式從雜亂的背景噪聲中有效提取其中的有用信息，大大減少了時頻圖原有的參數量，在一定程度上克服了傳統時頻分析中存在的缺陷。將該方法應用于大鼠嗅球LFP的研究，可提取與氣味刺激相關的事件特征，與支持向量機(support vector machine,SVM)算法相結合，初步實現了三種氣味刺激的分類。
1 分析方法
該方法的核心思想是用盡可能少的方程參數來表征原始時頻圖，以達到去除噪聲干擾，提取有用信息供進一步分析的目的。主要包括兩個步驟：① 構建神經信號基于小波變換的功率譜時頻圖；② 使用半橢球方程對時頻圖進行建模。
1.1 基于小波變換的功率譜時頻圖
小波變換是以某些特殊函數為基礎，將數據過程或數據系列變換為級數系列，以發現它的類似頻譜的特征，從而實現數據處理。這種方法是由法國工程師Morlet 在1974 年首先提出，經過之后數十年的研究，該方法建立了完整的數學化體系和扎實的理論基礎，并被廣泛應用于各種神經生理信號的分析[1-2, 10-13]。函數f(x)的連續小波變換可定義為
| ${{W}_{f}}\left( a,b \right)=\frac{1}{\sqrt{\left| a \right|}}\int_{-\infty }^{+\infty }{f\left( x \right)\bar{\Psi }}\left( \frac{x-b}{a} \right)dx,$ |
式中Ψ(x)是母小波，(x)是Ψ(x)的共軛形式，參數a 和b 分別被稱為尺度因子和平移因子。由式(1)可知，連續小波變換是空間(時間)和尺度的局部變換，因而能有效地從信號中提取信息。通過伸縮和平移等運算功能可對函數或信號進行多尺度的細化分析，解決了傅里葉變換不能解決的許多困難。
由連續小波變換得到的小波功率譜可定義為
| $TF\left( a,b \right)={{\left| {{W}_{f}}\left( a,b \right) \right|}^{2}}$ |
式中TF(a,b)是函數f(x)在尺度a和位移b時所具有的能量。嚴格地說，連續小波變換提供了時間尺度分析，而不是時間頻率分析。然而，使用適當的尺度頻率變換，依然可以得到一個近似于時間頻率分析的工具。與短時傅里葉變換相比，小波時頻圖的窗函數具有可變的特性。對于較高的頻率，時間窗會變小，使時間分析精度提高，而對于較低的頻率，時間窗會變大，保留了完整的低頻信息。
小波變換的關鍵是要選擇合適的母小波[1]。復Morlet小波是一種使用廣泛的基本小波,已被眾多研究證明，適用于電生理信號的時頻分析[1-3, 10-11, 13]。復Morlet小波w(t,f)對于某頻率的表達式為
| $w\left( t,f \right)=Aexp(2i\pi ft)exp({{t}^{2}}/2{{\sigma }_{i}}^{2})$ |
式中A為歸一化因子，其目的是使小波能量為1；f與σt(時域分辨率)是需要適當選擇的參數，兩者的乘積決定了小波的周期數，同一復Morlet小波家族具有恒定的f與σt乘積，本研究所使用的小波家族定義為2πfσt=7(見圖 1a)。復Morlet小波w(t,f)是一種復值調制的高斯(Gaussian)函數，在時域和頻域上都具有高斯分布，其傅里葉頻譜的標準差滿足
| ${{\sigma }_{f}}=\frac{1}{2\pi {{\sigma }_{t}}},$ |
特別的，在本研究中，
| ${{\sigma }_{f}}=\frac{1}{2\pi {{\sigma }_{t}}}=\frac{f}{7},$ |
因此，不同頻率f所對應的σt和σf不相同，即在整個時頻平面上具有可變的時頻分辨率：在高頻區能提供高的時間分辨率，在低頻區能提供高的頻率分辨率(見圖 1b)。
圖1
復Morlet小波的特性
(a)復Morlet小波的實部(公式3)，左：f/
(a) real part of equation 3,with
由于神經信號在低頻段具有的能量比高頻段大很多，為了在同一比例下觀察場電位信號在低頻段和高頻段能量隨時間的變化，由連續小波變換獲得的信號功率時頻圖需要標準化
| ${{z}_{ft}}=\frac{{{c}_{ft}}-{{m}_{f}}}{{{\delta }_{f}}},$ |
式中cft為f(t)在頻率f及時間t時的小波功率譜，mf和δf分別為刺激前1 s內信號在頻率f的小波功率譜的平均值和標準差。
圖 2為一段人造信號的原始波形及其基于復Morlet小波變換的時頻圖。該人造信號(見圖 2a)模擬了神經信號在各頻率段的成分組成，具有1個5 Hz的基礎振蕩，并在0.5~0.7 s和1.5 ~1.7 s范圍內分別具有一個β段(15~40 Hz)和σ段(40~90 Hz)的突變。圖 2(b)為該信號的原始時頻圖，由于低頻段基礎振蕩的能量大，使得高頻段的突變不明顯，經公式(6)歸一化之后，能量的突變點被凸顯出來，而恒定的能量波段被弱化。
圖2
人造信號及其小波時頻圖
(a)人造信號原始波形；(b)人造信號原始小波時頻圖；(c)經過歸一化之后的小波時頻圖
Figure2. An artificial signal and its wavelet time-frequency maps(a) the waveform of an artificial signal； (b) the raw wavelet time-frequency map； (c) the normalized wavelet time-frequency map
1.2 基于半橢球方程的時頻圖建模
原始時頻圖包含的數據量非常大，很難直接提取信號的時頻特征。在本研究中，采用了建模的方法進行信息提取和降緯處理。該方法的核心思想是利用事先預定的參數方程φ對原始時頻圖進行擬合，主要包含時頻圖窗口化、求和與擬合三個步驟。
首先將原始時頻圖分解成一系列重疊的窗口，以確定需要建模的區域。在經行窗口化時，窗口的尺寸必須滿足三點要求：① 由于原始時頻圖是通過連續小波變換獲得，在進行窗口化時必須使窗口的尺寸能夠自動適應小波變換的時頻分辨率，即L/σt=H/σf。其中L和H分別為窗口的長度和高度，σt和σf分別為小波變換的時域標準差和頻域標準差。 ② 要選擇合適的窗口寬度，在充分體現頻率特征的基礎上凸顯能量的突變，一般選擇的窗口需要包含當前頻率3~5 個振蕩周期。計算中，我們選取的值為4。 ③ 為了使位于不同位置的窗口具有相同的評判標準，窗口的面積必須相同。
因此，對于特定的頻率f，窗口的長度L=4/f，即一個窗口包含當前頻率的4個振蕩周期，則L/σt=8π/7，由此換算出窗口的高度H=8πf/49。用這個值去計算原始時頻圖上每個點，就能獲得一系列面積為3249π的重疊窗口。然后逐一計算每個窗口中的能量總和，窗口總能量值體現了該窗口內場電位的振蕩強度。對窗口wi，總能量值Si計算公式如式7所示，然后找出所有大于預設閾值的窗口。
| ${{S}_{i}}=\sum\limits_{t,f\in {{w}_{i}}}{{{z}_{ft}}}$ |
最后，在選中的窗口中使用參數方程進行擬合。本研究選擇半橢球方程φ作為擬合方程并定義為
| $\varphi \left( f,t \right)=\left\{ \begin{align} & a\sqrt{1-v},0 ＜v＜1 \\ & 0,v>1 \\ \end{align} \right.$ |
其中v=(f/μf)2/lf2+(t/μt)/lt2，μf和μt為該橢球在時頻圖上的位置，a為橢球的極半徑(幅度)，lf和lt為半橢球的赤道半徑。圖 3為一個半橢球方程的示意圖。該方程共有5個參數，并且需要滿足以下條件：
(1)μf ＞ 0，μt ＞ 0，且使半橢球位于所選窗口的中心；
(2)0 ＜ lt ＜ L,0 ＜ lf ＜ H。
圖3
由5個參數進行表征的半橢球方程示意圖
Figure3.
Half-ellipsoid model with its five parameters
在擬合的過程中，要使建模誤差函數C達到最小，函數C定義為
| $C=\frac{1}{2}\sum\limits_{t,f\in {{w}_{i}}}{({{z}_{ft}}-\varphi \left( f,t \right)),}$ |
圖 4為采自麻醉大鼠在香芹酮刺激下嗅球的局部場電位β段信號的小波時頻圖以及建模后的結果。從圖中我們可以發現，用5個半橢球模型就能將原始時頻圖的主要能量變化從大量參數和背景噪聲中提取出來。
圖4
嗅球場電位信號β段時頻圖和半橢球方程建模結果
左：歸一化后的小波時頻圖；中：半橢球建模的2D表示；右：半橢球建模的3D表示
Figure4. Time-frequency map (left) of β range of signal from rat’s olfactory bulb and its modeling results in two (middle) and three (right) dimensions2 大鼠嗅球場電位記錄和氣味識別分析
氣味分子的識別是由嗅黏膜上的一千多種氣味受體完成的，一種氣味受體只能與特定分子結構的氣味相結合[15-17]。表達同種氣味受體的感覺神經元投射到嗅球表面的一個或幾個嗅小球中，因此特定的氣味能夠使嗅球產生特定時空特性的活動[18-22]。研究者們發現，氣味刺激前后嗅球僧帽細胞層神經電活動十分強烈，LFP振蕩會明顯增加，其幅值和波形的變化受到氣味種類的調控 [9, 23-25]。
2.1 實驗設計和信號記錄
實驗對象SD大鼠購于浙江省醫學科學院，質量是200~250 g。大鼠進行電極植入手術流程是在腹腔注射10% 的水合氯醛(4 mL/kg)麻醉后，移除單側嗅球上覆蓋的頭骨，利用油壓微推進器(Narishige,Japan)將自制的16通道直徑為65 μm鎳鎘合金微絲(AM system,USA)電極陣列，植入到嗅球僧帽細胞層，平均深度約為200~300 μm。由于嗅球僧帽細胞放電幅值較大，能提示電極已植入到目標細胞層，故我們選取氣味刺激時能記錄到細胞鋒電位信號的通道，分析該通道的低頻場電位信號。
使用濃度為10 mmol/L的異丁醇、苯甲醚和香芹酮三種氣味作為刺激物，單次刺激給予體積為13 mL、時間為5 s的氣味，刺激間隔為60 s。多通道信號采集裝置為商用OmniPlex數據采集系統(Plexon,USA)，具有20倍的前置放大器，16位的A/D轉換器，能以最高頻率為40 kHz，5 000倍模擬放大，實現256個通道實時數據采集，并帶有用于低頻場電位和高頻鋒電位的頻帶隔離可編程硬件濾波器。信號采集裝置通過網線將采集到的數據發送到計算機，利用Matlab編譯環境實現數據的處理和分析。實驗的整體示意圖如圖 5所示。
圖5
實驗整體示意圖
Figure5.
Schematic of the experiment
2.2 特征提取和氣味分類
考慮到了小波變換和半橢球方程建模時可能造成的邊界效應，我們將記錄的嗅球場電位信號分成三個獨立的頻段分別進行變換。頻率的步進分別為0.25 Hz(θ段、β段)和1 Hz(γ段)。1.5~18.5 Hz：用于研究θ段的信號(5~12 Hz)；9~37.5 Hz：用于研究β段的信號(12~30 Hz)；25~85 Hz：用于研究γ段的信號(30~75 Hz)。經過半橢球方程建模之后，原始信號可以用一系列的半橢球的參數來表征。我們定義了以下參數作為分類的特征值：F1，橢球的個數；F2，橢球的最大幅值; F3，橢球的幅值總和。
因此，每個頻段能獲得3個特征值，每段信號能獲得9個特征值。這些特征值經過歸一化之后，用于后續的基于SVM算法的氣味分類。 本文采用了SVM算法[26]對氣味信號行進分類識別。三種氣味共218組時頻圖向量(苯甲醚73組，異丁醇72組，香芹酮73組)，隨機選取一半作為樣本A，另一半作為樣本B，進行分類測試，最終識別結果如表 1所示。二重交叉檢驗得到的總體準確率為79.4%，其中苯甲醚、異丁醇和香芹酮的檢測準確率分別為82.2%、80.6%和75.3%，其準確度和特性性均高于在以往研究中基于短時傅里葉變換的分析和分類方法[9]。
表1
氣味分類結果
Table1.
Results of odor discrimination
3 結論
本文介紹了一種新型的低頻神經信號的處理方法。該方法將基于小波變換的時頻圖分析方法與半橢球建模方法有機結合起來，能夠快速有效地提取神經信號能量變化的有用信息，也可用于比較大量數據在相同與不同刺激下的相同與不同之處，克服了傳統時頻分析方法背景噪聲大、參數多的缺陷。該方法目前存在的主要問題是建模算法比較復雜，當待處理數據段增長時，對計算機硬件配置要求較高，數據處理時間大幅增加，有待進一步的優化。
將該方法用于麻醉大鼠嗅球場電位的分析中，可用于提取不同氣味刺激下，嗅球的電生理響應特征，通過模式識別的方法，初步實現不同氣味刺激的分類，其準確率和特異性均高于基于短時傅里葉變換的分析和分類方法。
