對誘發電位(EP)信號中具有強脈沖過程的腦電圖(EEG)噪聲, 可以用α穩定分布模型來描述。基于分數低階矩對傳統的Cohen類時頻分布進行了改進, 得到了新的分數低階空間時頻分布(FLO-STFM), 據此提出了一種新的可在α穩定分布環境下工作的分數低階空間時頻欠定盲分離算法(FLO-TF-UBSS)。將該盲分離算法應用到EP信號的提取, 仿真實驗結果表明所提出的盲分離算法能較好地在EEG噪聲環境下實現對EP信號的盲提取, 相關系數以及盲提取效果都優于基于二階的TF-UBSS算法。
引用本文： 龍俊波, 汪海濱, 查代奉. 基于分數低階空間時頻矩陣的腦電誘發電位盲提取. 生物醫學工程學雜志, 2015, 32(2): 269-274. doi: 10.7507/1001-5515.20150049 復制
引言
誘發電位(evoked potential, EP)是大腦中樞神經系統對外界因素刺激所產生的一種條件反應，通常可以作為腦病分析的重要依據。在對EP信號的提取過程中，往往不可避免地會混入一種可以被認為是噪聲的自發腦電圖(electroencephalograph, EEG)信號，由于EEG噪聲有很強的脈沖性，所以一般性噪比非常小。盲源分離(blind source separation, BSS)作為一種重要的手段已經被廣泛應用到EP信號的提取中，主要提取方法有二階盲識別算法[1]、旋轉不變性信號參數估計算法、多未知信號提取算法、矩陣聯合近似對角化算法、快速獨立分量分析法[2]和擴展信息最大化算法[3]。文獻[1]研究了各種算法對軀體感覺誘發電位的提取，并進行了綜合比較，確定了二階盲識別算法優勢。Ahmadian等[4]提出了一種受限盲提取方法，實現對預備潛在電位的提取，體現了一定的性能優勢。文獻[3]則利用自回歸最小均方算法和基于小波的主分量分析提出了一種新算法，實現了對眼部活動的偽量消除。
在對EP信號提取建模的過程中，一些學者認為含有EEG噪聲或者其它噪聲的EP信號不可以簡單地用高斯模型來描述，加速度撞擊、缺氧等特殊試驗時EEG噪聲通常脈沖性很強，具有重拖尾及無限方差性特征，其概率密度函數服從α穩定分布[5-8]。大量數據表明其特征指數為1.06~1.94，如果采用傳統的二階高斯模型來描述，算法性能將大大減退甚至可能失效。因此，需要選擇具有更好適用性的新的α穩定分布統計模型。最近幾年，國外針對α穩定分布模型的EP信號提取算法已經有一定的研究，如文獻[5]提出了一種基于分數低階最小分散系數的盲分離算法。國內相關研究也取得了一定的成果。如文獻[6]提出了基于BOREL譜測度的ICA分析盲分離算法；林政劍等[7]提出一種基于分數低階矩共變的盲分離算法；文獻[8]在最小均方誤差算法的基礎上提出了一種數據復用自適應最小平均P范數算法。這些算法都局限于時域數據分析。最近，基于時頻域的盲分離算法開始受到了關注。Aissa-El-Bey等[9]提出了基于短時傅里葉變換(short-time Fourier transform, STFT)時頻和偽魏格納-維納時頻的欠定盲分離算法(time-frequency underdetermined blind source separation, TF-UBSS)，實現了時頻域非相交源的盲分離。基于時頻分布的欠定混疊盲分離算法在文獻[10]中被提出。文獻[9-10]算法都是基于二階空間時頻分布矩陣(spatial time-frequency matrix，STFD)提出的，當信號或噪聲為α穩定分布過程時，顯然算法會退化或者失效。因此，本文基于分數低階矩，對魏格納-維納時頻分布(Wigner-Ville distribution，WVD)、偽魏格納-維納時頻分布(pseudo WVD，PWVD)進行了改進，得到了適合α穩定分布過程的分數低階魏格納-維納時頻分布(fractional lower order WVD，FLO-WVD)和分數低階偽魏格納-維納時頻分布(fractional lower order PWVD，FLO-PWVD)時頻分布，并以此提出了分數低階空間時頻分布矩陣(fractional lower order spatial time-frequency matrix, FLO-STFM)的概念，結合TF-UBSS算法得到了一種新的可應用到穩定分布環境的分數低階時頻欠定盲分離算法(fractional lower order TF-UBSS, FLO-TF-UBSS)算法。最后，總結了算法的步驟，分析了算法的性能，并將其實際應用到穩定分布EEG噪聲下的EP盲提取。
1 分數低階時頻盲提取算法
1.1 算法模型及時頻分布
含EEG噪聲的EP信號采集通常用腦電采集儀獲取，多通道傳感器同時采集到信號和噪聲，信號和噪聲線性混合后可以用觀測數據模型表示為：
| ${y_t}=As\left(t \right)=\sum\limits_{i=1}^N {{a_i}{s_i}\left(t \right)} $ |
我們假設接收通道傳感器的數目為M，信號和噪聲的個數為N，在處理過程中我們把噪聲也當做信號來對待，假設信號和噪聲相互獨立。那么，觀測信號y(t)=[y1(t), y2(t), …, yM(t)]T，線性混合矩陣A=[a1, a2, …, aN]，aN=[aN1, aN2, …, aNM]T, 所以A是一個M×N的矩陣，N個信號和噪聲s(t)=[s1(t), ss(t), …, sN(t)]T，假設N個信號和噪聲在時頻域不相交。基于這些假設，從信號和噪聲的時頻域獲得信號和噪聲各自的時頻分布，再利用時頻重構還原信號，實現對EP信號的提取。
信號x(t)的WVD及PWVD可以表示為：
| $WV{D_x}\left({t, f} \right)=\int\limits_{-\infty }^\infty {x\left({t + \tau /2} \right){x^*}\left({t-\tau /2} \right){{\rm{e}}^{-j2\pi fr}}{\rm{d}}\tau } $ |
| $WV{D_x}\left({t, f} \right)=\int\limits_{-\infty }^\infty {h\left(\tau \right)x\left({t + \tau /2} \right){x^*}\left({t-\tau /2} \right){{\rm{e}}^{-j2\pi fr}}{\rm{d}}\tau } $ |
由于分數低階過程沒有有限的二階矩，那么式(2)、(3)時頻分布定義失去了意義，為此我們根據分數低階矩的定義x〈P〉=｜x｜psign(x)進行改進。其中當x＞0時sign(x)=1，當x=0時sign(x)=0，當x＜0時sign(x)=-1，x-〈P〉=(x*)〈P〉=(x〈P〉)*，*表示共軛，p＜α/2，α是穩定分布信號或噪聲的特征指數。改進后的得到的分數低階WVD和PWVD計算式分別表示為：
| $FLOPWV{D_x}\left({t, f} \right)=\int\limits_{-\infty }^\infty {{x^{\left\langle P \right\rangle }}\left({t + \tau /2} \right){x^{-\left\langle P \right\rangle }}\left({t-\tau /2} \right){{\rm{e}}^{-j2\pi fr}}{\rm{d}}\tau } $ |
| $FLOPWV{D_x}\left({t, f} \right)=\int\limits_{-\infty }^\infty {h\left(\tau \right)\left({t + \tau /2} \right){x^{-\left\langle P \right\rangle }}\left({t-\tau /2} \right){{\rm{e}}^{-j2\pi fr}}{\rm{d}}\tau } $ |
1.2 分數低階時頻盲提取算法原理
觀測數據y(t)=[y1(t), y2(t), …, yM(t)]T的FLO-STFM可以定義為：
| ${\Gamma _Y}\left({t, f} \right)=\left[{\begin{array}{*{20}{c}} {FLOPWV{D_{{y_1}{y_2}\left({t, f} \right)}}} & \cdots & {FLOPWV{D_{{y_1}{y_M}}}\left({t, f} \right)}\ \vdots & \ddots & \vdots \ {FLOPWV{D_{{y_M}{y_1}}}\left({t, f} \right)} & \cdots & {FLOPWVD{y_M}{y_M}\left({t, f} \right)} \end{array}} \right]=A{\Gamma _S}\left({t, f} \right){A^H}$ |
其中FLOPWVDyiyj(t, f)是觀測陣元yi(t)與觀測陣元yj(t)之間的FLO-PWVD, ΓS(t, f)是信號和噪聲的FLO-PWVD，ΓYY(t, f)為混合FLO-STFM。
首先，考慮到小噪聲的影響，我們可以設置一個閾值ε(通常取0.05)對FLO-STFM中所有的時頻分布點(tp, fp)進行篩選，便可以得到屬于ΓY(t, f)的時頻支撐域，計算式為：
| $\frac{{{\Gamma _Y}{{\left({{t_p}, {f_p}} \right)}^{\left\langle P \right\rangle }}}}{{{{\max }_f}\left\{ {{\Gamma _Y}{{\left({{t_p}, {f_p}} \right)}^{\left\langle P \right\rangle }}} \right\}}} > \varepsilon $ |
對每一個(tp, fp)進行搜索，如果滿足式(7)條件，則認為屬于ΓY(t, f)的時頻支撐域，并保留該點(tp, fp)的值，如果不滿足式(7)條件，則令該點(tp, fp)值為零，這個過程可以稱為噪聲“軟”閾值濾波過程。
經過式(7)閾值后，我們再對ΓY(t, f)進行白化處理，由于二階信號有有限的方差，通常可以利用協方差矩陣的求逆方根來得到白化矩陣，針對含有分數低級EEG噪聲的EP信號，其方差無限，所以，白化方法不可以采用傳統的方法。y(t)=[y1(t), y2(t), …, yM(t)]T的歸一化協方差矩陣定義為：
| $\Phi=\frac{{y\left(t \right){y^T}\left(t \right)}}{{M \cdot {\rm{trace}}\left\{ {y\left(t \right){y^T}\left(t \right)/M} \right\}}}$ |
式(8)中trace為矩陣的跡，對Φ進行奇異值分解，我們可以得到Φ=UΛ2UT, 那么預白化矩陣W=Λ-1UT, 令白化后的FLO-STFM為, 那么有：
| ${\underline \Gamma _{YY}}\left({{t_p}, {f_p}} \right)=W{\Gamma _{YY}}\left({{t_p}, {f_p}} \right){W^H}$ |
PWVD存在一定的交叉項，即EEG噪聲和EP信號之間的時頻點以及EP信號與EP信號之間的互時頻點。為了消除這些交叉項的影響，我們可以利用式(10)準則進行處理。
| $\frac{{{\rm{trace}}\left\{ {{{\underline \Gamma }_{YY}}\left({{t_p}, {f_p}} \right)} \right\}}}{{{{\underline \Gamma }_{YY}}{{\left({{t_p}, {f_p}} \right)}^{\left\langle P \right\rangle }}}} < {\varepsilon _0}$ |
ε0為一設點的門限值(一般取0.85)，如果滿足式(10)條件，那么認為(tp, fp)為互時頻點，否則認為是自時頻點。經過式(10)篩選后，DYY(tp, fp)就只剩下那些含有EP信號和EEG噪聲的自時頻點(tp, fp)，由于每一個信號源都有一定的時頻聚集性，信號源之間又不相交，如兩個信號源的時頻分布中，信號1時頻集中在Ω1，信號2時頻為Ω2，那么，所有的信號源可以表示為，對每一個信號源的時頻聚集區域Ωi有：
| ${\Gamma _{YY}}\left({{t_p}, {f_p}} \right)={\Gamma _{{S_i}{S_i}}}\left({{t_p}, {f_p}} \right){a_i}a_i^H, \forall \left({{t_p}, {f_p}} \right) \in {\Omega _i}$ |
那么，我們只要得到EP信號以及EEG噪聲的向量矩陣a(tp, fp)，便可以還原出信號和噪聲。文獻[10]采用了張量分解的方法對混合矩陣進行估計，而文獻[9]選擇了在已知信號源的數目情況下，采用K階聚集算法進行估計，對每一個(tp, fp)∈Ωi時頻點，我們根據文獻[9]所提出的算法進行改進，空間向量a(tp, fp)的計算式為：
| $a\left({{t_p}, {f_p}} \right)=\frac{{{\rm{diag}}\left\{ {{{\underline \Gamma }_{YY}}\left({{t_p}, {f_p}} \right)} \right\}}}{{{\rm{diag}}{{\left\{ {{{\underline \Gamma }_{YY}}\left({{t_p}, {f_p}} \right)} \right\}}^{\left\langle P \right\rangle }}}}$ |
得到一系列的{a(tp, fp)｜(tp, fp)∈Ω}后，因為具有相同主特征向量的所有信號自項屬于同一個信號源，具有相同的空間方向，所以對這些向量利用文獻[9]中的聚類算法可以進行分類，就可以得到每一個源信號的向量{ai(tp, fp)｜(tp, fp)∈Ωi}，代入式(11)就可以得到對應的ΓYY(tp, fp)，把所有屬于該信號源的時頻點ΓYY(tp, fp)集合后便可以得到該信號源單獨的時頻分布估計, 計算式為：
| ${\mathop {FLOPWVD}\limits^ \wedge _{{S_i}}}\left({t, f} \right)={\rm{trace}}\left\{ {{{\underline \Gamma }_{YY}}\left({t, f} \right)} \right\}, \left({t, f} \right) \in {\Omega _i}$ |
最后，用式(14)對得到的進行時頻重構，便可以得到對應的源信號：
| $s\left(t \right)=\frac{1}{{{s^*}\left(0 \right)}}\int\limits_{-\infty }^\infty {{{\mathop {FLOPWVD}\limits^ \wedge }_S}\left({\frac{t}{2}, f} \right){{\rm{e}}^{j2\pi ft}}{\rm{d}}f} $ |
1.3 分數低階時頻盲提取算法步驟
(1)利用式(5)和式(6)計算觀測信號矩陣的分數低階空間時頻矩陣ΓYY(t, f)。
(2)先利用式(7)對ΓYY(t, f)進行噪聲閾值，再基于歸一化協方差矩陣利用式(8)和式(9)對信號進行預白化處理，得到白化后的矩陣。
(3)利用式(10)對EP信號或EEG噪聲的自時頻點進行選擇，消除交叉項的干擾。
(4)基于K階聚集算法式(12)計算每一個源信號的空間向量{ai(tp, fp)｜(tp, fp)∈Ωi}，并利用式(13)計算每一個信號的單獨分數低階時頻分布。
(5)基于得到的利用式(14)進行時頻重構，還原出EP信號或EEG噪聲si(t)。
1.4 評價指標
為了評定所提出的分數低階空間時頻盲提取方法和已有的基于二階的盲提取算法分離效果，我們選擇評價指標為相關系數，計算式如式(15)所示，其中si(n)為源信號，為盲分離估計得到的信號。
| $\rho \left({{s_i}, {{\hat s}_j}} \right)=\frac{{\left| {\sum\limits_{n=1}^N {{s_i}\left(n \right){{\hat s}_j}\left(n \right)} } \right|}}{{\sqrt {\sum\limits_{n=1}^N {s_i^2\left(n \right)\sum\limits_{n=1}^N {\hat s_j^2\left(n \right)} } } }}$ |
由于不能使用信噪比，所以我們使用廣義信噪比(generalized signal-to-noise ratio，GSNR)為EP信號和EEG噪聲(SαS)的信噪比：
| ${\rm{GSNR}}=10\log \left(E \right)\left\{ {{{\left| {s\left(t \right)} \right|}^2}} \right\}/\left. {{\gamma ^a}} \right)$ |
2 計算機仿真
誘發電位信號通常以加速度撞擊或缺氧窒息試驗條件誘發并通過腦電圖儀采集獲得，取500個點的EP信號(采樣頻率1 000 Hz)，SαS過程的EEG噪聲用計算機模擬產生，如圖 1所示，用一個3×2的隨機矩陣A(如23式)把EP信號和EEG噪聲進行線性混合，同時加入20 dB的高斯隨機噪聲，混合后的觀測信號如圖 2所示。
| $A=\left[{\begin{array}{*{20}{c}} {2.8396} & {0.8914}\ {-2.0080} & {-0.2330}\ {-1.7048} & {-0.5537} \end{array}} \right]$ |
圖1
EP信號和EEG噪聲
Figure1.
EP signal and EEG noise
圖2
混合后的觀測信號
Figure2.
Waveform of the observation point
分別采用TF-UBSS算法和本文所提出的FLO-TF-UBSS算法進行仿真試驗，具體實驗如下。
2.1 EP信號和EEG噪聲盲提取
本實驗取EEG噪聲特征指數α=1.5，EP信號和EEG噪聲GSNR=30 dB，分別用TF-UBSS和FLO-TF-UBSS(取P=1.1)對EP信號和EEG噪聲進行盲提取，得到的EP信號和EEG噪聲信號如圖 3、圖 4所示。
圖3
TF-UBSS算法盲提取效果
Figure3.
Blind extraction effect of TF-UBSS algorithm
圖4
FLO-TF-UBSS算法盲提取效果
Figure4.
Blind extraction effect of FLO-TF-UBSS algorithm
2.2 不同GSNR下相關系數比較
本實驗取EEG噪聲特征指數α=1.5，取EP信號和EEG噪聲GSNR從10~30 dB進行混合，然后分別用TF-UBSS和FLO-TF-UBSS(P=1.1)進行盲提取，計算得到的相關系數如表 1所示。
表1
不同的GSNR下相關系數比較
Table1.
Comparison of correlation coefficient between different GSNR
2.3 不同α下相關系數比較
本實驗取GSNR=30 dB，EEG噪聲的特征指數α從1.06~1.94變化，對混合后的觀測信號使用TF-UBSS和FLO-TF-UBSS(P=1.1)算法對EP信號和EEG噪聲進行分離，用相關系數衡量分離效果，仿真效果如圖 5、圖 6所示。
圖5
不同α下EP相關系數
Figure5.
EP correlation coefficient on different α
圖6
不同α下EEG相關系數
Figure6.
EEG correlation coefficient on different α
從圖 3、圖 4兩種算法的提取效果來看，由于SαS穩定分布EEG噪聲過程具有較強的脈沖性，導致TF-UBSS提取算法(圖 3)性能明顯退化，而圖 4中FLO-TF-UBSS算法卻體現了優良的性能，分離出來的EP信號及EEG噪聲與圖 1中的原信號幾乎相同，從效果觀察FLO-TF-UBSS算法優于TF-UBSS算法。另外，從表 1中可以看出，不管是低信噪比還是高信噪比情況下，所提出的FLO-TF-UBSS算法相關系數都接近1，且相關系數始終都大于TF-UBSS算法的相關系數，而TF-UBSS算法相關系數則偏差較大，充分體現了FLO-TF-UBSS盲提取效果良好。從圖 5、圖 6中可以看出，不同的特征指數條件下，FLO-TF-UBSS算法提取后的EP和EEG相關系數都近似為1，而TF-UBSS算法隨著特征指數從1.06~1.94逐漸接近1，但始終小于FLO-TF-UBSS提取算法的相關系數。所以，本文所提出的FLO-TF-UBSS算法在對EP的提取上，不論是從分離效果還是從穩定度來看都優于已有的TF-UBSS算法。
3 結論
針對特殊激勵條件，如速度撞擊或缺氧環境，EEG信號具有強脈沖、重拖尾及無限方差性。因此，α穩定分布可以更好地描述這些特殊情況下EP信號中存在的EEG信號。本研究改進了基于二階的傳統時頻分布，得到了基于分數低階矩的時頻分布，提出了FLO-STFM的概念，并基于FLO-STFM提出了一種適合α穩定分布環境的FLO-TF-UBSS盲提取算法。仿真結果表明，FLO-TF-UBSS盲提取算法能有效地分離出EEG信號，相比已有的TF-UBSS盲提取算法有更好的分離效果，更接近原信號。不同GSNR下相關系數比較分析比較表明，FLO-TF-UBSS算法相關系數始終接近1，而TF-UBSS算法相關系數則有一定的偏差。不同的特征指數α下，FLO-TF-UBSS算法的EP和EEG相關系數約等于1，而TF-UBSS算法隨著特征指數變化而變化。總之，所提出的FLO-TF-UBSS算法優于TF-UBSS算法，是一種具有良好韌性的盲提取方法。
引言
誘發電位(evoked potential, EP)是大腦中樞神經系統對外界因素刺激所產生的一種條件反應，通常可以作為腦病分析的重要依據。在對EP信號的提取過程中，往往不可避免地會混入一種可以被認為是噪聲的自發腦電圖(electroencephalograph, EEG)信號，由于EEG噪聲有很強的脈沖性，所以一般性噪比非常小。盲源分離(blind source separation, BSS)作為一種重要的手段已經被廣泛應用到EP信號的提取中，主要提取方法有二階盲識別算法[1]、旋轉不變性信號參數估計算法、多未知信號提取算法、矩陣聯合近似對角化算法、快速獨立分量分析法[2]和擴展信息最大化算法[3]。文獻[1]研究了各種算法對軀體感覺誘發電位的提取，并進行了綜合比較，確定了二階盲識別算法優勢。Ahmadian等[4]提出了一種受限盲提取方法，實現對預備潛在電位的提取，體現了一定的性能優勢。文獻[3]則利用自回歸最小均方算法和基于小波的主分量分析提出了一種新算法，實現了對眼部活動的偽量消除。
在對EP信號提取建模的過程中，一些學者認為含有EEG噪聲或者其它噪聲的EP信號不可以簡單地用高斯模型來描述，加速度撞擊、缺氧等特殊試驗時EEG噪聲通常脈沖性很強，具有重拖尾及無限方差性特征，其概率密度函數服從α穩定分布[5-8]。大量數據表明其特征指數為1.06~1.94，如果采用傳統的二階高斯模型來描述，算法性能將大大減退甚至可能失效。因此，需要選擇具有更好適用性的新的α穩定分布統計模型。最近幾年，國外針對α穩定分布模型的EP信號提取算法已經有一定的研究，如文獻[5]提出了一種基于分數低階最小分散系數的盲分離算法。國內相關研究也取得了一定的成果。如文獻[6]提出了基于BOREL譜測度的ICA分析盲分離算法；林政劍等[7]提出一種基于分數低階矩共變的盲分離算法；文獻[8]在最小均方誤差算法的基礎上提出了一種數據復用自適應最小平均P范數算法。這些算法都局限于時域數據分析。最近，基于時頻域的盲分離算法開始受到了關注。Aissa-El-Bey等[9]提出了基于短時傅里葉變換(short-time Fourier transform, STFT)時頻和偽魏格納-維納時頻的欠定盲分離算法(time-frequency underdetermined blind source separation, TF-UBSS)，實現了時頻域非相交源的盲分離。基于時頻分布的欠定混疊盲分離算法在文獻[10]中被提出。文獻[9-10]算法都是基于二階空間時頻分布矩陣(spatial time-frequency matrix，STFD)提出的，當信號或噪聲為α穩定分布過程時，顯然算法會退化或者失效。因此，本文基于分數低階矩，對魏格納-維納時頻分布(Wigner-Ville distribution，WVD)、偽魏格納-維納時頻分布(pseudo WVD，PWVD)進行了改進，得到了適合α穩定分布過程的分數低階魏格納-維納時頻分布(fractional lower order WVD，FLO-WVD)和分數低階偽魏格納-維納時頻分布(fractional lower order PWVD，FLO-PWVD)時頻分布，并以此提出了分數低階空間時頻分布矩陣(fractional lower order spatial time-frequency matrix, FLO-STFM)的概念，結合TF-UBSS算法得到了一種新的可應用到穩定分布環境的分數低階時頻欠定盲分離算法(fractional lower order TF-UBSS, FLO-TF-UBSS)算法。最后，總結了算法的步驟，分析了算法的性能，并將其實際應用到穩定分布EEG噪聲下的EP盲提取。
1 分數低階時頻盲提取算法
1.1 算法模型及時頻分布
含EEG噪聲的EP信號采集通常用腦電采集儀獲取，多通道傳感器同時采集到信號和噪聲，信號和噪聲線性混合后可以用觀測數據模型表示為：
| ${y_t}=As\left(t \right)=\sum\limits_{i=1}^N {{a_i}{s_i}\left(t \right)} $ |
我們假設接收通道傳感器的數目為M，信號和噪聲的個數為N，在處理過程中我們把噪聲也當做信號來對待，假設信號和噪聲相互獨立。那么，觀測信號y(t)=[y1(t), y2(t), …, yM(t)]T，線性混合矩陣A=[a1, a2, …, aN]，aN=[aN1, aN2, …, aNM]T, 所以A是一個M×N的矩陣，N個信號和噪聲s(t)=[s1(t), ss(t), …, sN(t)]T，假設N個信號和噪聲在時頻域不相交。基于這些假設，從信號和噪聲的時頻域獲得信號和噪聲各自的時頻分布，再利用時頻重構還原信號，實現對EP信號的提取。
信號x(t)的WVD及PWVD可以表示為：
| $WV{D_x}\left({t, f} \right)=\int\limits_{-\infty }^\infty {x\left({t + \tau /2} \right){x^*}\left({t-\tau /2} \right){{\rm{e}}^{-j2\pi fr}}{\rm{d}}\tau } $ |
| $WV{D_x}\left({t, f} \right)=\int\limits_{-\infty }^\infty {h\left(\tau \right)x\left({t + \tau /2} \right){x^*}\left({t-\tau /2} \right){{\rm{e}}^{-j2\pi fr}}{\rm{d}}\tau } $ |
由于分數低階過程沒有有限的二階矩，那么式(2)、(3)時頻分布定義失去了意義，為此我們根據分數低階矩的定義x〈P〉=｜x｜psign(x)進行改進。其中當x＞0時sign(x)=1，當x=0時sign(x)=0，當x＜0時sign(x)=-1，x-〈P〉=(x*)〈P〉=(x〈P〉)*，*表示共軛，p＜α/2，α是穩定分布信號或噪聲的特征指數。改進后的得到的分數低階WVD和PWVD計算式分別表示為：
| $FLOPWV{D_x}\left({t, f} \right)=\int\limits_{-\infty }^\infty {{x^{\left\langle P \right\rangle }}\left({t + \tau /2} \right){x^{-\left\langle P \right\rangle }}\left({t-\tau /2} \right){{\rm{e}}^{-j2\pi fr}}{\rm{d}}\tau } $ |
| $FLOPWV{D_x}\left({t, f} \right)=\int\limits_{-\infty }^\infty {h\left(\tau \right)\left({t + \tau /2} \right){x^{-\left\langle P \right\rangle }}\left({t-\tau /2} \right){{\rm{e}}^{-j2\pi fr}}{\rm{d}}\tau } $ |
1.2 分數低階時頻盲提取算法原理
觀測數據y(t)=[y1(t), y2(t), …, yM(t)]T的FLO-STFM可以定義為：
| ${\Gamma _Y}\left({t, f} \right)=\left[{\begin{array}{*{20}{c}} {FLOPWV{D_{{y_1}{y_2}\left({t, f} \right)}}} & \cdots & {FLOPWV{D_{{y_1}{y_M}}}\left({t, f} \right)}\ \vdots & \ddots & \vdots \ {FLOPWV{D_{{y_M}{y_1}}}\left({t, f} \right)} & \cdots & {FLOPWVD{y_M}{y_M}\left({t, f} \right)} \end{array}} \right]=A{\Gamma _S}\left({t, f} \right){A^H}$ |
其中FLOPWVDyiyj(t, f)是觀測陣元yi(t)與觀測陣元yj(t)之間的FLO-PWVD, ΓS(t, f)是信號和噪聲的FLO-PWVD，ΓYY(t, f)為混合FLO-STFM。
首先，考慮到小噪聲的影響，我們可以設置一個閾值ε(通常取0.05)對FLO-STFM中所有的時頻分布點(tp, fp)進行篩選，便可以得到屬于ΓY(t, f)的時頻支撐域，計算式為：
| $\frac{{{\Gamma _Y}{{\left({{t_p}, {f_p}} \right)}^{\left\langle P \right\rangle }}}}{{{{\max }_f}\left\{ {{\Gamma _Y}{{\left({{t_p}, {f_p}} \right)}^{\left\langle P \right\rangle }}} \right\}}} > \varepsilon $ |
對每一個(tp, fp)進行搜索，如果滿足式(7)條件，則認為屬于ΓY(t, f)的時頻支撐域，并保留該點(tp, fp)的值，如果不滿足式(7)條件，則令該點(tp, fp)值為零，這個過程可以稱為噪聲“軟”閾值濾波過程。
經過式(7)閾值后，我們再對ΓY(t, f)進行白化處理，由于二階信號有有限的方差，通常可以利用協方差矩陣的求逆方根來得到白化矩陣，針對含有分數低級EEG噪聲的EP信號，其方差無限，所以，白化方法不可以采用傳統的方法。y(t)=[y1(t), y2(t), …, yM(t)]T的歸一化協方差矩陣定義為：
| $\Phi=\frac{{y\left(t \right){y^T}\left(t \right)}}{{M \cdot {\rm{trace}}\left\{ {y\left(t \right){y^T}\left(t \right)/M} \right\}}}$ |
式(8)中trace為矩陣的跡，對Φ進行奇異值分解，我們可以得到Φ=UΛ2UT, 那么預白化矩陣W=Λ-1UT, 令白化后的FLO-STFM為, 那么有：
| ${\underline \Gamma _{YY}}\left({{t_p}, {f_p}} \right)=W{\Gamma _{YY}}\left({{t_p}, {f_p}} \right){W^H}$ |
PWVD存在一定的交叉項，即EEG噪聲和EP信號之間的時頻點以及EP信號與EP信號之間的互時頻點。為了消除這些交叉項的影響，我們可以利用式(10)準則進行處理。
| $\frac{{{\rm{trace}}\left\{ {{{\underline \Gamma }_{YY}}\left({{t_p}, {f_p}} \right)} \right\}}}{{{{\underline \Gamma }_{YY}}{{\left({{t_p}, {f_p}} \right)}^{\left\langle P \right\rangle }}}} < {\varepsilon _0}$ |
ε0為一設點的門限值(一般取0.85)，如果滿足式(10)條件，那么認為(tp, fp)為互時頻點，否則認為是自時頻點。經過式(10)篩選后，DYY(tp, fp)就只剩下那些含有EP信號和EEG噪聲的自時頻點(tp, fp)，由于每一個信號源都有一定的時頻聚集性，信號源之間又不相交，如兩個信號源的時頻分布中，信號1時頻集中在Ω1，信號2時頻為Ω2，那么，所有的信號源可以表示為，對每一個信號源的時頻聚集區域Ωi有：
| ${\Gamma _{YY}}\left({{t_p}, {f_p}} \right)={\Gamma _{{S_i}{S_i}}}\left({{t_p}, {f_p}} \right){a_i}a_i^H, \forall \left({{t_p}, {f_p}} \right) \in {\Omega _i}$ |
那么，我們只要得到EP信號以及EEG噪聲的向量矩陣a(tp, fp)，便可以還原出信號和噪聲。文獻[10]采用了張量分解的方法對混合矩陣進行估計，而文獻[9]選擇了在已知信號源的數目情況下，采用K階聚集算法進行估計，對每一個(tp, fp)∈Ωi時頻點，我們根據文獻[9]所提出的算法進行改進，空間向量a(tp, fp)的計算式為：
| $a\left({{t_p}, {f_p}} \right)=\frac{{{\rm{diag}}\left\{ {{{\underline \Gamma }_{YY}}\left({{t_p}, {f_p}} \right)} \right\}}}{{{\rm{diag}}{{\left\{ {{{\underline \Gamma }_{YY}}\left({{t_p}, {f_p}} \right)} \right\}}^{\left\langle P \right\rangle }}}}$ |
得到一系列的{a(tp, fp)｜(tp, fp)∈Ω}后，因為具有相同主特征向量的所有信號自項屬于同一個信號源，具有相同的空間方向，所以對這些向量利用文獻[9]中的聚類算法可以進行分類，就可以得到每一個源信號的向量{ai(tp, fp)｜(tp, fp)∈Ωi}，代入式(11)就可以得到對應的ΓYY(tp, fp)，把所有屬于該信號源的時頻點ΓYY(tp, fp)集合后便可以得到該信號源單獨的時頻分布估計, 計算式為：
| ${\mathop {FLOPWVD}\limits^ \wedge _{{S_i}}}\left({t, f} \right)={\rm{trace}}\left\{ {{{\underline \Gamma }_{YY}}\left({t, f} \right)} \right\}, \left({t, f} \right) \in {\Omega _i}$ |
最后，用式(14)對得到的進行時頻重構，便可以得到對應的源信號：
| $s\left(t \right)=\frac{1}{{{s^*}\left(0 \right)}}\int\limits_{-\infty }^\infty {{{\mathop {FLOPWVD}\limits^ \wedge }_S}\left({\frac{t}{2}, f} \right){{\rm{e}}^{j2\pi ft}}{\rm{d}}f} $ |
1.3 分數低階時頻盲提取算法步驟
(1)利用式(5)和式(6)計算觀測信號矩陣的分數低階空間時頻矩陣ΓYY(t, f)。
(2)先利用式(7)對ΓYY(t, f)進行噪聲閾值，再基于歸一化協方差矩陣利用式(8)和式(9)對信號進行預白化處理，得到白化后的矩陣。
(3)利用式(10)對EP信號或EEG噪聲的自時頻點進行選擇，消除交叉項的干擾。
(4)基于K階聚集算法式(12)計算每一個源信號的空間向量{ai(tp, fp)｜(tp, fp)∈Ωi}，并利用式(13)計算每一個信號的單獨分數低階時頻分布。
(5)基于得到的利用式(14)進行時頻重構，還原出EP信號或EEG噪聲si(t)。
1.4 評價指標
為了評定所提出的分數低階空間時頻盲提取方法和已有的基于二階的盲提取算法分離效果，我們選擇評價指標為相關系數，計算式如式(15)所示，其中si(n)為源信號，為盲分離估計得到的信號。
| $\rho \left({{s_i}, {{\hat s}_j}} \right)=\frac{{\left| {\sum\limits_{n=1}^N {{s_i}\left(n \right){{\hat s}_j}\left(n \right)} } \right|}}{{\sqrt {\sum\limits_{n=1}^N {s_i^2\left(n \right)\sum\limits_{n=1}^N {\hat s_j^2\left(n \right)} } } }}$ |
由于不能使用信噪比，所以我們使用廣義信噪比(generalized signal-to-noise ratio，GSNR)為EP信號和EEG噪聲(SαS)的信噪比：
| ${\rm{GSNR}}=10\log \left(E \right)\left\{ {{{\left| {s\left(t \right)} \right|}^2}} \right\}/\left. {{\gamma ^a}} \right)$ |
2 計算機仿真
誘發電位信號通常以加速度撞擊或缺氧窒息試驗條件誘發并通過腦電圖儀采集獲得，取500個點的EP信號(采樣頻率1 000 Hz)，SαS過程的EEG噪聲用計算機模擬產生，如圖 1所示，用一個3×2的隨機矩陣A(如23式)把EP信號和EEG噪聲進行線性混合，同時加入20 dB的高斯隨機噪聲，混合后的觀測信號如圖 2所示。
| $A=\left[{\begin{array}{*{20}{c}} {2.8396} & {0.8914}\ {-2.0080} & {-0.2330}\ {-1.7048} & {-0.5537} \end{array}} \right]$ |
圖1
EP信號和EEG噪聲
Figure1.
EP signal and EEG noise
圖2
混合后的觀測信號
Figure2.
Waveform of the observation point
分別采用TF-UBSS算法和本文所提出的FLO-TF-UBSS算法進行仿真試驗，具體實驗如下。
2.1 EP信號和EEG噪聲盲提取
本實驗取EEG噪聲特征指數α=1.5，EP信號和EEG噪聲GSNR=30 dB，分別用TF-UBSS和FLO-TF-UBSS(取P=1.1)對EP信號和EEG噪聲進行盲提取，得到的EP信號和EEG噪聲信號如圖 3、圖 4所示。
圖3
TF-UBSS算法盲提取效果
Figure3.
Blind extraction effect of TF-UBSS algorithm
圖4
FLO-TF-UBSS算法盲提取效果
Figure4.
Blind extraction effect of FLO-TF-UBSS algorithm
2.2 不同GSNR下相關系數比較
本實驗取EEG噪聲特征指數α=1.5，取EP信號和EEG噪聲GSNR從10~30 dB進行混合，然后分別用TF-UBSS和FLO-TF-UBSS(P=1.1)進行盲提取，計算得到的相關系數如表 1所示。
表1
不同的GSNR下相關系數比較
Table1.
Comparison of correlation coefficient between different GSNR
2.3 不同α下相關系數比較
本實驗取GSNR=30 dB，EEG噪聲的特征指數α從1.06~1.94變化，對混合后的觀測信號使用TF-UBSS和FLO-TF-UBSS(P=1.1)算法對EP信號和EEG噪聲進行分離，用相關系數衡量分離效果，仿真效果如圖 5、圖 6所示。
圖5
不同α下EP相關系數
Figure5.
EP correlation coefficient on different α
圖6
不同α下EEG相關系數
Figure6.
EEG correlation coefficient on different α
從圖 3、圖 4兩種算法的提取效果來看，由于SαS穩定分布EEG噪聲過程具有較強的脈沖性，導致TF-UBSS提取算法(圖 3)性能明顯退化，而圖 4中FLO-TF-UBSS算法卻體現了優良的性能，分離出來的EP信號及EEG噪聲與圖 1中的原信號幾乎相同，從效果觀察FLO-TF-UBSS算法優于TF-UBSS算法。另外，從表 1中可以看出，不管是低信噪比還是高信噪比情況下，所提出的FLO-TF-UBSS算法相關系數都接近1，且相關系數始終都大于TF-UBSS算法的相關系數，而TF-UBSS算法相關系數則偏差較大，充分體現了FLO-TF-UBSS盲提取效果良好。從圖 5、圖 6中可以看出，不同的特征指數條件下，FLO-TF-UBSS算法提取后的EP和EEG相關系數都近似為1，而TF-UBSS算法隨著特征指數從1.06~1.94逐漸接近1，但始終小于FLO-TF-UBSS提取算法的相關系數。所以，本文所提出的FLO-TF-UBSS算法在對EP的提取上，不論是從分離效果還是從穩定度來看都優于已有的TF-UBSS算法。
3 結論
針對特殊激勵條件，如速度撞擊或缺氧環境，EEG信號具有強脈沖、重拖尾及無限方差性。因此，α穩定分布可以更好地描述這些特殊情況下EP信號中存在的EEG信號。本研究改進了基于二階的傳統時頻分布，得到了基于分數低階矩的時頻分布，提出了FLO-STFM的概念，并基于FLO-STFM提出了一種適合α穩定分布環境的FLO-TF-UBSS盲提取算法。仿真結果表明，FLO-TF-UBSS盲提取算法能有效地分離出EEG信號，相比已有的TF-UBSS盲提取算法有更好的分離效果，更接近原信號。不同GSNR下相關系數比較分析比較表明，FLO-TF-UBSS算法相關系數始終接近1，而TF-UBSS算法相關系數則有一定的偏差。不同的特征指數α下，FLO-TF-UBSS算法的EP和EEG相關系數約等于1，而TF-UBSS算法隨著特征指數變化而變化。總之，所提出的FLO-TF-UBSS算法優于TF-UBSS算法，是一種具有良好韌性的盲提取方法。
