基于克里金插值的EEG二維成像最優(yōu)化分析

蔡 靖， 胡紀(jì)鋒， 劉鋒華， 蒙堅(jiān)發(fā)

（吉林大學(xué)儀器科學(xué)與電氣工程學(xué)院，長(zhǎng)春 130061）

0 引 言

腦電信號(hào)（Electroencephalography EEG）是腦神經(jīng)細(xì)胞群產(chǎn)生的微弱的非平穩(wěn)偽隨機(jī)生物電信號(hào)，含有豐富的大腦活動(dòng)狀態(tài)信息。隨著近年信息技術(shù)的發(fā)展，運(yùn)用腦電信號(hào)進(jìn)行腦電成像分析并應(yīng)用于臨床研究越來(lái)越成為研究的重點(diǎn)，腦電成像可以直觀顯示大腦活動(dòng)狀態(tài)，對(duì)于腦血管疾病、腦積水、腦腫瘤等病灶的定位［1-3］以及癲癇、中風(fēng)等疾病發(fā)作時(shí)腦電成像特點(diǎn)［4-5］的臨床醫(yī)學(xué)研究有很重要的意義。

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

在以上研究的基礎(chǔ)上，本文對(duì)基于高斯模型變異函數(shù)普通克里金插值算法進(jìn)行優(yōu)化分析，在運(yùn)用數(shù)據(jù)擬合理論的基礎(chǔ)上建立了變異函數(shù)各參量之間的函數(shù)關(guān)系，并結(jié)合模型的均方誤差給出了優(yōu)化的普通克里金變異函數(shù)模型取值規(guī)范。

1 腦電成像算法分析

1.1 普通克里金插值

克里金插值是Matheron首先提出［10-12］的一種對(duì)有限區(qū)域內(nèi)部的區(qū)域化變量進(jìn)行無(wú)偏最優(yōu)估計(jì)的空間插值算法。

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

式中：x為位置坐標(biāo)；h為距離；Z（x＋h）、Z（x）分別為x＋h與x位置處的實(shí)測(cè)值；γ（h）為距離h對(duì)應(yīng)的半變異函數(shù)。

對(duì)于需要插值獲得的單個(gè)待測(cè)樣本點(diǎn)x0，由n個(gè)樣本點(diǎn)的實(shí)測(cè)值Z（xi）可得到待測(cè)點(diǎn)Z*（x0）的估計(jì)值：

式中，λ為權(quán)重系數(shù)。則估計(jì)值的方差為：

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

對(duì)于多個(gè)待測(cè)點(diǎn)，通常采用估計(jì)值方差的平均值來(lái)衡量插值效果。估計(jì)值方差的平均值為：

式中：N為待測(cè)點(diǎn)數(shù)；n為已知實(shí)際測(cè)量點(diǎn)的個(gè)數(shù)。估計(jì)值方差的平均值作為腦電插值成像精度的評(píng)判標(biāo)準(zhǔn)，當(dāng)估計(jì)值方差的平均值＝0時(shí)即認(rèn)為是最優(yōu)的腦電成像。

1.2 半變異函數(shù)模型與優(yōu)化關(guān)系確定

式（5）中，半變異函數(shù)γ（h）表征了插值區(qū)域的空間變異結(jié)構(gòu)，一個(gè)插值區(qū)域的半變異函數(shù)定義為：

實(shí)際運(yùn)用中為了簡(jiǎn)化計(jì)算，通常選用下列半變異函數(shù)模型。

高斯模型：

指數(shù)模型：

球型模型：

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

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

的C0與ɑ的值即為高斯模型的最優(yōu)參數(shù)。

1.3 相關(guān)尺度a與塊金值最優(yōu)值獲取

對(duì)于n個(gè)樣本點(diǎn)，以任意兩個(gè)樣本點(diǎn)之間的距離h作為橫坐標(biāo)，由式（6）得出的C（h）值作為縱坐標(biāo)，進(jìn)行二次曲線擬合，即可得到取樣區(qū)域的近似半變異函數(shù)，擬合結(jié)果如圖1所示。

圖1 二次擬合曲線

擬合多項(xiàng)式如下：

聯(lián)立求解式（7）、（11）式可得ɑ與C0的關(guān)系：

對(duì)于不同的距離h，ɑ和C0有相似的關(guān)系，取所有h的平均值h0代入上式，得到最優(yōu)的ɑ和C0的關(guān)系（見圖2）。

圖2 ɑ-C0關(guān)系曲線

聯(lián)立式（10）、（12）得：

即可求出C0與ɑ的最優(yōu)值。

由式（13）得出的ɑ與C0參數(shù)進(jìn)行優(yōu)化普通克里金插值，并用此插值結(jié)果進(jìn)行最優(yōu)化腦電成像。

2 腦電信號(hào)采集系統(tǒng)

如圖3所示為腦電信號(hào)采集系統(tǒng)結(jié)構(gòu)圖，主要由采集電極、前置濾波放大、模數(shù)轉(zhuǎn)換、微處理器、藍(lán)牙通信和上位機(jī)信號(hào)處理等單元組成。

圖3 系統(tǒng)結(jié)構(gòu)圖

采集電極安放位置參考國(guó)際10-20系統(tǒng)電極放置法，8個(gè)采集電極Fp1、Fp2、Fz、Cz、P3、P4、O1、O2同時(shí)采集腦電信號(hào)，信號(hào)經(jīng)由前置濾波放大單元濾去高頻段的噪聲，并將信號(hào)放大。模數(shù)轉(zhuǎn)換單元將8通道的模擬腦電信號(hào)轉(zhuǎn)換為數(shù)字信號(hào)，微處理器單元負(fù)責(zé)接收此數(shù)字信號(hào)并將其通過藍(lán)牙通信傳輸?shù)缴衔粰C(jī)中進(jìn)行信號(hào)的處理。

3 實(shí)驗(yàn)及結(jié)果分析

3.1 原始腦電信號(hào)

由采集系統(tǒng)采集Fp1、Fp2、Fz、Cz、P3、P4、O1、O2點(diǎn)的相對(duì)于耳垂的原始電位值腦電信號(hào)如圖4所示。

3.2 測(cè)量電極坐標(biāo)

為進(jìn)行腦電2D成像，以國(guó)際10-20電極安放標(biāo)準(zhǔn)中的CZ為極點(diǎn)，以CZ點(diǎn)到眉間的弧線為極軸建立極坐標(biāo)系，極徑為電極安放點(diǎn)到原點(diǎn)的大腦弧面距離，極角為電極安放點(diǎn)與極點(diǎn)相連線段在極點(diǎn)所在平面的投影與極軸的夾角。

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

建立平面直角坐標(biāo)系，將極坐標(biāo)點(diǎn)轉(zhuǎn)換為直角坐標(biāo)點(diǎn)并平移到第1象限內(nèi)。即可得到電極安放位置坐標(biāo)如圖5所示。

圖5 測(cè)量電極坐標(biāo)圖

3.3 半變異函數(shù)模型選擇

普通克里金插值中變異函數(shù)模型的選取對(duì)插值的無(wú)偏性具有重要意義。圖6～8分別為變異函數(shù)模型為高斯模型、指數(shù)模型、球形模型時(shí)腦電成像結(jié)果及對(duì)應(yīng)的誤差分布圖。模型待測(cè)點(diǎn)估計(jì)值方差分布情況如圖9所示，不同半變異函數(shù)模型插值點(diǎn)誤差曲線如圖10所示。

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

圖7 指數(shù)模型腦電成像及誤差圖

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

由圖可見，高斯模型每個(gè)待測(cè)點(diǎn)的估計(jì)值誤差均在零附近，且平均誤差明顯低于其他2種模型，因此本文進(jìn)行普通克里金插值優(yōu)化分析時(shí)選用的半變異函數(shù)模型為高斯模型。

3.4 最優(yōu)值確定

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

圖9 方差散點(diǎn)圖

圖10 不同半變異函數(shù)模型插值點(diǎn)誤差曲線

圖11 平均方差

由式（13）可得參數(shù)最優(yōu)解。圖12中曲面與曲線交點(diǎn)即反映了當(dāng)平均誤差趨近于0時(shí)的相關(guān)尺度ɑ與塊金值C0的最優(yōu)值。

圖12 交點(diǎn)圖

3.5 結(jié)果分析

圖13為在正常靜默狀態(tài)下，相關(guān)尺度ɑ與塊金值C0取位于最優(yōu)解附近的值時(shí)測(cè)試者的2D腦電成像情況及對(duì)應(yīng)的誤差分布圖。

圖13 腦電成像及誤差圖

均方差值如表1所示。

以上實(shí)驗(yàn)的基礎(chǔ)上采用逐次逼近的方法，可得均方差趨近于0時(shí)的最優(yōu)的ɑ和C0值。

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

在均方差趨近于零的約束條件下，對(duì)高斯模型半變異函數(shù)的參數(shù)ɑ和C0的函數(shù)關(guān)系進(jìn)行最優(yōu)解求取，實(shí)現(xiàn)了基于高斯模型半變異函數(shù)的普通克里金插值的優(yōu)化分析。實(shí)驗(yàn)結(jié)果表明，采用逐步逼近的方法平均方差趨近于零。

表1 均方差值表

圖14 最優(yōu)腦電成像及誤差圖

4 結(jié) 語(yǔ)

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