基于克里金插值的EEG二維成像最優化分析

蔡 靖， 胡紀鋒， 劉鋒華， 蒙堅發

（吉林大學儀器科學與電氣工程學院，長春 130061）

0 引 言

腦電信號（Electroencephalography EEG）是腦神經細胞群產生的微弱的非平穩偽隨機生物電信號，含有豐富的大腦活動狀態信息。隨著近年信息技術的發展，運用腦電信號進行腦電成像分析并應用于臨床研究越來越成為研究的重點，腦電成像可以直觀顯示大腦活動狀態，對于腦血管疾病、腦積水、腦腫瘤等病灶的定位［1-3］以及癲癇、中風等疾病發作時腦電成像特點［4-5］的臨床醫學研究有很重要的意義。

Lee等［6］開發了一臺帶有IBM PC AT的電腦腦電圖成像系統，采用反距離加權插值的方法，使每個插值點的權重與其最近的4個點的距離成反比，實現了腦電成像。Yang等［7］借助混合效應多項式回歸模型的統計分析方法，成功實現了用高密度腦電圖來描述次最大肌肉收縮時腦電成像中相應動態源的變化。Paul等［8］對于癲癇病癥進行了研究，結果表明，在大多數情況下使用k-最近鄰點插值法可以有效提高靈敏度，并且誤差也較小。Taran等［9］使用IMF1的徑向基核函數對不同運動圖像任務的腦電信號進行分類，該方法與現有的方法相比，在分類準確度，靈敏性等方面展現了更好的性能。

在以上研究的基礎上，本文對基于高斯模型變異函數普通克里金插值算法進行優化分析，在運用數據擬合理論的基礎上建立了變異函數各參量之間的函數關系，并結合模型的均方誤差給出了優化的普通克里金變異函數模型取值規范。

1 腦電成像算法分析

1.1 普通克里金插值

克里金插值是Matheron首先提出［10-12］的一種對有限區域內部的區域化變量進行無偏最優估計的空間插值算法。

克里金插值法中的普通克里金法具有下列約束條件：

式中：x為位置坐標；h為距離；Z（x＋h）、Z（x）分別為x＋h與x位置處的實測值；γ（h）為距離h對應的半變異函數。

對于需要插值獲得的單個待測樣本點x0，由n個樣本點的實測值Z（xi）可得到待測點Z*（x0）的估計值：

式中，λ為權重系數。則估計值的方差為：

式中：hi為x0到xi的距離；L為拉格朗日乘數。

對于多個待測點，通常采用估計值方差的平均值來衡量插值效果。估計值方差的平均值為：

式中：N為待測點數；n為已知實際測量點的個數。估計值方差的平均值作為腦電插值成像精度的評判標準，當估計值方差的平均值＝0時即認為是最優的腦電成像。

1.2 半變異函數模型與優化關系確定

式（5）中，半變異函數γ（h）表征了插值區域的空間變異結構，一個插值區域的半變異函數定義為：

實際運用中為了簡化計算，通常選用下列半變異函數模型。

高斯模型：

指數模型：

球型模型：

式中：C0為塊金值，是由樣本誤差和短距離的變異性引起的半變異函數的偏差；ɑ為相關尺度，表示當樣本之間的距離大于等于此距離時，各樣本相互獨立。C0＋C1＝var［Z（x）］，稱為基臺值。

可以得出C0與ɑ的取值決定了方差均值的大小。經過比較，本文選用高斯模型，在高斯模型下使

的C0與ɑ的值即為高斯模型的最優參數。

1.3 相關尺度a與塊金值最優值獲取

對于n個樣本點，以任意兩個樣本點之間的距離h作為橫坐標，由式（6）得出的C（h）值作為縱坐標，進行二次曲線擬合，即可得到取樣區域的近似半變異函數，擬合結果如圖1所示。

圖1 二次擬合曲線

擬合多項式如下：

聯立求解式（7）、（11）式可得ɑ與C0的關系：

對于不同的距離h，ɑ和C0有相似的關系，取所有h的平均值h0代入上式，得到最優的ɑ和C0的關系（見圖2）。

圖2 ɑ-C0關系曲線

聯立式（10）、（12）得：

即可求出C0與ɑ的最優值。

由式（13）得出的ɑ與C0參數進行優化普通克里金插值，并用此插值結果進行最優化腦電成像。

2 腦電信號采集系統

如圖3所示為腦電信號采集系統結構圖，主要由采集電極、前置濾波放大、模數轉換、微處理器、藍牙通信和上位機信號處理等單元組成。

圖3 系統結構圖

采集電極安放位置參考國際10-20系統電極放置法，8個采集電極Fp1、Fp2、Fz、Cz、P3、P4、O1、O2同時采集腦電信號，信號經由前置濾波放大單元濾去高頻段的噪聲，并將信號放大。模數轉換單元將8通道的模擬腦電信號轉換為數字信號，微處理器單元負責接收此數字信號并將其通過藍牙通信傳輸到上位機中進行信號的處理。

3 實驗及結果分析

3.1 原始腦電信號

由采集系統采集Fp1、Fp2、Fz、Cz、P3、P4、O1、O2點的相對于耳垂的原始電位值腦電信號如圖4所示。

3.2 測量電極坐標

為進行腦電2D成像，以國際10-20電極安放標準中的CZ為極點，以CZ點到眉間的弧線為極軸建立極坐標系，極徑為電極安放點到原點的大腦弧面距離，極角為電極安放點與極點相連線段在極點所在平面的投影與極軸的夾角。

將人的大腦近似為一個球面，假設某兩個電極Fp1與Fz的球面坐標分別為r1、θ1、?1與r1、θ2、?2，那么兩者之間的球面距離，即弧長l可以用如下方程組求解：

建立平面直角坐標系，將極坐標點轉換為直角坐標點并平移到第1象限內。即可得到電極安放位置坐標如圖5所示。

圖5 測量電極坐標圖

3.3 半變異函數模型選擇

普通克里金插值中變異函數模型的選取對插值的無偏性具有重要意義。圖6～8分別為變異函數模型為高斯模型、指數模型、球形模型時腦電成像結果及對應的誤差分布圖。模型待測點估計值方差分布情況如圖9所示，不同半變異函數模型插值點誤差曲線如圖10所示。

圖6 高斯模型腦電成像及誤差圖

圖7 指數模型腦電成像及誤差圖

圖8 球形模型腦電成像及誤差圖

由圖可見，高斯模型每個待測點的估計值誤差均在零附近，且平均誤差明顯低于其他2種模型，因此本文進行普通克里金插值優化分析時選用的半變異函數模型為高斯模型。

3.4 最優值確定

當相關尺度ɑ與塊金值C0取值不同時，由式（4）可得平均方差變化曲線如圖11所示。

圖9 方差散點圖

圖10 不同半變異函數模型插值點誤差曲線

圖11 平均方差

由式（13）可得參數最優解。圖12中曲面與曲線交點即反映了當平均誤差趨近于0時的相關尺度ɑ與塊金值C0的最優值。

圖12 交點圖

3.5 結果分析

圖13為在正常靜默狀態下，相關尺度ɑ與塊金值C0取位于最優解附近的值時測試者的2D腦電成像情況及對應的誤差分布圖。

圖13 腦電成像及誤差圖

均方差值如表1所示。

以上實驗的基礎上采用逐次逼近的方法，可得均方差趨近于0時的最優的ɑ和C0值。

圖14為身體健康的成年男性測試者在ɑ和C0取最優值即均方差趨近于0時所繪制的最優2D腦電成像與誤差分布圖：

在均方差趨近于零的約束條件下，對高斯模型半變異函數的參數ɑ和C0的函數關系進行最優解求取，實現了基于高斯模型半變異函數的普通克里金插值的優化分析。實驗結果表明，采用逐步逼近的方法平均方差趨近于零。

表1 均方差值表

圖14 最優腦電成像及誤差圖

4 結 語

本文提出了一種基于克里金插值的EEG二維成像最優化分析方法。通過仿真比較，使用高斯半變異函數模型的EEG二維成像的方差均值最小。通過擬合高斯半變異函數中的協方差函數，得到相關尺度ɑ與塊金值C0的函數關系。方差均值是誤差評估標準，當方差均值接近零時，ɑ與C0將取得最優值。通過將ɑ與C0的最優值應用于EEG二維成像，實現了最優EEG二維成像。該方法提高了EEG二維成像的準確性，使得分析EEG信號并從EEG中獲取生理和疾病信息變得更加容易。