2020研究生資料建模C題總結-腦電訊號分析判別
前言
2020年研究生數學建模大賽上月結束，最終的比賽成績還沒有公佈，現在是2020年10月，在此總結一下建模大賽的選題與成果。記錄過去，展望未來。
我們團隊選擇C題，題目內容如下：
腦電訊號按其產生的方式可分為誘發腦電訊號和自發腦電訊號。誘發腦電訊號是通過某種外界刺激使大腦產生電位變化從而形成的腦電活動；自發腦電訊號是指在沒有外界特殊刺激下，大腦自發產生的腦電活動。
- P300電訊號：是誘發腦電訊號的一種，在小概率刺激發生後300毫秒範圍左右出現的一個正向波峰，如下圖所示。
它不受刺激物理特性影響，與知覺或認知心理活動有關，與注意、記憶、智慧等加工過程密切相關。被測者觀察字元矩陣(圖1-a)，並記下左上角的字元，然後按順序閃爍字元矩陣行列(圖1-b)的同時，記錄被測者的P300腦電訊號。- 睡眠腦電：是自發腦電訊號的一種，能夠反映身體狀態的自身變化，也是用來診斷和治療相關疾病的重要依據。根據睡眠狀態進一步分為清醒期、睡眠I期、睡眠II期、深度睡眠和快速眼動期，腦電訊號如圖2所示。
P300腦電資料由5位被測者的20次視覺刺激實驗的腦電資料組成，每次實驗的字元迴圈閃爍5輪，記錄頭頂20個位置的腦電訊號，腦電訊號記錄位置如下圖所示。
睡眠腦電資料3000各睡眠腦電樣本，由狀態標籤、Alpha、Beta、Theta和Delta波在“8-13Hz”，“14-25Hz”，“4-7Hz”和“0.5-4Hz”頻率範圍內的能量佔比。
本文主要解決以下3個問題：
1、對P300腦電資料去除冗餘資訊降維。
2、對不同視覺刺激下的腦電波進行分類，判斷閃爍的行或列中是否包含目標字元，最終確定目標字元。
資料預處理
建立資料字典
P300腦機介面資料包含訓練用的12個已知字元的腦電資料，和測試用的8個未知字元的腦電資料。將訓練資料和測試資料，以字元為鍵腦電資料為值，組成字典，如表1所示。
| 字典 | “Sx_train_data_dict” | “Sx_test_data_dict” | “Sx_train_event_dict” | “Sx_test_event_dict” |
|---|---|---|---|---|
| 鍵 | 目標字元 | 待測字元 | 目標字元 | 待測字元 |
| 值 | 字元對應的P300資料表 | 字元對應的P300資料表 | 字元對應的事件標籤表 | 字元對應的事件標籤表 |
生成資料索引表
計算索引
每個資料表首行首列為字元序號num
r o w = n u m / / 6 + 1 row = num // 6 +1 row=num//6+1
c o l = n u m % 7 + 1 col = num \% 7 +1 col=num%7+1
建立索引表
針對訓練集和測試集的不同字元對應的事件標籤表，去除識別符號和結束標籤所在的行，生成與P300資料表對應的索引表，STrait(t)表示一次閃爍刺激記錄的腦電資料：
第1列：字元矩陣某次STrait(t)發生位置，記錄行號或列號。
第2列：字元矩陣某次STrait(t)開始時對應P300資料表的位置索引
第3列：字元矩陣某次STrait(t)結束時對應P300資料表的位置索引
第4列：字元矩陣某次STrait(t)記錄的樣本數
第5列：1（閃爍中包含目標字元）或0（閃爍中不包含目標字元）
資料規整
由於每次STrait(t)的P300資料取樣情況不同，導致每次取樣資料量存在波動。對所有STrait(t)的資料取樣量進行統計分析，取樣量最小值為3，最大值為46，平均數為39，下四分位數為41，中位數為41，上四分位數為43。因此，刪除取樣資料量小於40的索引記錄。
資料切片
索引表中每行資料記錄了每個STrait(t)的閃爍的行/列、開始位置索引、結束位置索引、取樣樣本量和包含目標字元標誌。逐行遍歷索引表，將之對應的P300資料表切片，獲得每次STrait(t)的P300資料。
問題一：資料去冗餘
每個被試者的腦電訊號，對於每個目標字元閃爍事件記錄的腦電訊號共20個通道，大約3000多條資料。原始腦電資料量較大，資料已經在正負類比例上產生較大差距，正例較少導致訓練結果一般。另一方面，原資料包含的資訊過多，部分通道間的資訊有較強相關性，造成資料冗餘，對分類產生一定的干擾。
度量指標
皮爾遜相關係數
皮爾遜相關係數是最早由統計學家卡爾·皮爾遜設計的統計指標，是研究變數之間線性相關程度的量，一般用字母 r 表示。
相關係數是按積差方法計算，同樣以兩變數與各自平均值的離差為基礎，通過兩個離差相乘來反映兩變數之間相關程度；著重研究線性的單相關係數。依據相關現象之間的不同特徵，其統計指標的名稱有所不同，將反映兩變數間線性相關關係的統計指標稱為相關係數。
兩變數間的皮爾遜相關係數可通過以下公式計算：
Cov(X,Y)為X與Y的協方差，Var[X]為X的方差，Var[Y]為Y方差。對於變數X和Y，當相關係數為0時，X和Y兩變數無關係；當X的值增大（減小），Y值增大（減小），兩個變數為正相關，相關係數在0.00與1.00之間；當X的值增大（減小），Y值減小（增大），兩個變數為負相關，相關係數在-1.00與0.00之間。相關係數的絕對值越大，相關性越強，相關係數越接近於1或-1，相關度越強，相關係數越接近於0，相關度越弱。通常情況下通過以下取值範圍判斷變數的相關強度：0.8-1.0 為極強相關；0.6-0.8為強相關；0.4-0.6為中等程度相關；0.2-0.4為弱相關；0.0-0.2為極弱相關或無相關。
皮爾遜相關係數實現過程
對於每個表格迴圈讀取每行腦電資料，並儲存為最終二維陣列。接著將每個通道作為一組向量，求取兩兩向量間的皮爾遜相關係數，生成大小為[20,20]的相關係數矩陣，可將其存為表格以便觀察，如圖5所示，是第5位受試者訓練資料相關係數矩陣。行列索引表示通道識別符號，為該通道的資料向量，矩陣具有沿對角線的對稱性，所以只觀察下三角區即可。圖示中的資料說明，不同通道向量之間的相關性不同，有正相關有負相關，有相關性較強的，也有無關的。相關係數越高，說明這兩個通道的腦波的變化特徵越相似，相關性極強的可只用一個通道表徵這種特徵，另一個通道是冗餘的。
主成分分析
主成分分析（Principal Component Analysis，PCA），是一種統計方法。通過正交變換將一組可能存在相關性的變數轉換為一組線性不相關的綜合變數，轉換後的這組變數叫主成分，主成分數量較原變數儘量少。在研究多變數任務或者做資料分析處理時，變數個數太多會增加任務複雜性，希望變數個數較少而得到的資訊較多。圖6為主成分分析原理示意圖。
主成分分析實現過程
對於腦電訓練資料的陣列，與每個通道的均值做差，然後計算協方差矩陣，並求解協方差矩陣的特徵值和特徵向量。將特徵值從大到小進行排序，作為一個特徵陣列儲存，以便之後使用。按照順序，越靠前特徵值越大，包含的資訊量就越多，越靠後則資訊量越少。
結果融合
對於每位受試者，遍歷訓練資料中所有目標字元的子表，對其對應陣列計算相關係數矩陣和PCA的特徵排序陣列。對於相關係數矩陣，每位受試者設定不同的閾值，係數大於此值的，將其行列索引即通道識別符號進行儲存。對於每個子表，對篩選出的相關通道對，比較其在特徵排序陣列中的前後位置，若比較靠後的通道其索引>10，則說明此通道包含的資訊較少，且其資訊少於另一通道，則將此通道刪除，另一通道保留，依次迴圈所有通道對。
上述過程，得到特定受試者的12個子表的12個通道組合，將所有組合合成一維陣列，計算1-20通道識別符號出現的次數，若出現次數>10，則啟用此通道，最後得到篩選出的個人最優通道組合。
實驗中輸出的相關資料如下，圖7為受試者1的12個子表PCA特徵排序結果，可以看出從13通道以後排列各不相同，說明不同目標字元實驗中的腦波也有一定的變化，與實際情況相符。根據受試者生理特徵的獨特性，相關係數的選擇閾值設定也不同。閾值越高，說明對資料相關性或負相關性的要求越高，達到相關性要求的通道對就越少，自然從特徵排序中刪除的通道數目減少，反之在組合中留下的通道越多。但是由於選取通道個數有一定限制，我們將各受試者組合通道數控制在14和16之間，因此對於不同受試者選定了不同閾值。或者，根據全域性最優通道數的要求，可以調整。
所有被試者通道融合
在得出各受試者的最有通道組合後，採取同理的融合方式，將5個通道組合數組合成一維陣列，對通道識別符號計數，若大於3次，則啟用此通道。基於上述結果，全域性最優通道組合識別符號為下圖所示，共包含13個通道，對應的通道名稱為Fz、F3、F4、Cz、C3、C4、T7、T8、CP3、CP4、CP5、Pz、P3。有了此結果，可以通過給腦電資料降維，刪除冗餘通道，留下對應的主特徵，進一步提高分類準確率。
問題二：P300腦電訊號分類
實驗流程
以經過去冗餘處理後選取得13個通道資料作為分類輸入，通過多種模型得效能對比，選取其中最優的模型，實驗流程如下圖所示。
資料去噪——濾波
腦電訊號時非平穩得，對分類視覺刺激的準確率有很大影響。為提高模型識別準確率與穩定性，我們考慮對腦電訊號進行濾波處理。經驗模態分解(empirical mode decomposition，EMD) 演算法[1]是一種將隨機訊號分解成多個內涵模態分量（IMF），將非平穩訊號轉化為平穩訊號的演算法。在該問題中，我們對給出的非平穩腦電訊號序列，依次找出其極大值點與極小值點，通過樣條插值，擬合出上下包絡曲線，根據上下包絡曲線求出平均包絡曲線，將原訊號資料減去該平均包絡曲線便得出分解的第一個IMF分量，將上述過程進行N次迭代，得到N個IMF分量，直至擬合曲線無負區域性極小值與正區域性極大值，停止迭代，公式如下所示：
其中，x(t)代表原非穩定隨機訊號，Ci (t)表示經過第i次EMD分解之後的IMF分量，rn(t)代表原訊號被分解成一系列IMF分量後的剩餘部分。
在解決問題1時，我們對附件1中的訊號序列進行EMD分解，結果如下圖所示，其中imf1~7為從原始序列中分離出來的頻率依次從高至低的區域性特徵訊號，res為經過7次imf分量分離後的剩餘訊號分量。為減少高低頻噪聲對後續模型訓練的影響，我們將原始訊號分別剔除高低頻分量作為後續的訓練資料。
模型效果比較
輸入預處理完畢的Train，Test資料集。並將其劃分為10等分，訓練集佔2成，輪次為4。本文采用共計6個模型，進行分類識別比較，實驗如下所示：
通過實驗可以發現，具有較好精度效能的模型為Random Forest、Decision TreeGini、fastDTW-KNN。其中Decision TreeGini計算時間最短，但測試集精度最低；SVM演算法訓練集精度良好，但測試集精度較低且耗時較長；fastDTW-KNN效能相對穩定。
問題三：睡眠腦電資料分類
提供3000個睡眠腦電特徵樣本及其標籤資料，第一列為已知標籤，第二至五列為從原始時序列中計算得到的特徵引數，依次包括“Alpha”，“Beta”，“Theta”，“Delta”，分別對應了腦電訊號在“8-13Hz”，“14-25Hz”，“4-7Hz”和“0.5-4Hz”頻率範圍內的能量佔比，特徵引數單位為百分比。其中，清醒期共記錄633條，快速眼動期共599條，睡眠期（I期）共記錄512條，睡眠期（II期）共記錄604條，深睡眠期共記錄602條。
本問題旨在使用盡可能少的訓練樣本資料的基礎上實現相對較高的預測準確率，為此，我們將訓練集資料佔比從0.7、0.5、0.3、0.1，依次輸入模型訓練並測試，最終確定資料佔比和模型選擇。
模型效能對比
候選實驗模型在不同資料集比例下的準確率如下所示。
通過對比實驗可知，在不斷減少訓練集佔比的情況下，兩候選模型普遍出現過擬合現象。通過模型對驗證集和測試集分類的準確率結果中可知，隨機森林模型相對於SVM具有更加穩定的預測效能。
總結與感想
本次比賽，在現有的機器學習分類任務研究的基礎上開展了深入細緻的研究，通過P300腦電波資料，比較了經典機器學習分類演算法的優劣，以皮爾遜係數與PCA的融合策略進一步緩解了採集通道的資料冗餘問題，優化了KNN演算法以提升分類效率。
一次比賽，在這麼短時間內完成較為複雜的任務挺不容易的，感謝一路走來問我們提供幫助的同學朋友，特別要感謝我可愛的女朋友提供的幫助與鼓勵。雖不知最終結果如何，但這次經歷加深了我對機器學習演算法的理解和資料探勘工作的實踐，希望可以在資料探勘領域繼續走下去。