楊淳沨 吳國成 伍家松 姜龍玉 孔佑勇 舒華忠

(東南大學計算機網(wǎng)絡和信息集成教育部重點實驗室， 南京 210096)(東南大學計算機科學與工程學院， 南京 210096)(東南大學中法生物醫(yī)學信息研究中心, 南京 210096)

近年來，通過分析和處理大腦電信號來研究大腦連通性已成為腦工作原理和腦部重大疾病發(fā)生機制研究領域的熱點之一[1-2].WGCI (Wiener-Granger causality index)算法[3]、DTF (directed transfer function)算法[4]、PDC (partial directed coherence)算法[5]等都可有效應用于研究大腦功能連通性和效應連通性[6].這些算法都是基于隨機過程的自回歸(autoregressive, AR)模型的算法，模型參數(shù)(階數(shù)和系數(shù))的估計起著舉足輕重的作用.傳統(tǒng)的AR模型階數(shù)估計算法包括赤池信息準則(Akaike information criterion, AIC)[7]、最后預測誤差法(final prediction error, FPE)和最小描述長度法(minimum description length, MDL)[8]等.這些算法均假設不同信號具有相同階數(shù).然而，在實際情況中，不同信號各自的階數(shù)往往各不相同，因此這些算法會使模型項數(shù)過估計，從而導致后續(xù)的系數(shù)估計運算量加大且精準度不高.學者們繼而開發(fā)了更加有效的AR模型參數(shù)估計算法[9-11].其中，最佳參數(shù)搜索(optimal parameter search, OPS)算法[12]是一種魯棒性較高的算法，即使在存在噪聲且模型初始階數(shù)選取不當?shù)那闆r下，OPS算法仍能篩選出正確的模型項.但是，在干擾項較多的情況下，OPS算法的性能下降較快，從而會產(chǎn)生較大誤差，導致難以篩選出所有正確的模型項.

本文提出了一種基于gAIC (generalized AIC)[13]和滑動窗的自回歸模型參數(shù)估計算法.其中，模型的階數(shù)估計采用gAIC算法，與AIC相比，它能更有效地減少干擾項數(shù)目.對原始信號進行加窗處理，可以消除部分信號的不穩(wěn)定性.實驗結(jié)果表明，相比傳統(tǒng)算法，所提算法在模型參數(shù)估計方面具有更好的魯棒性.

1 模型參數(shù)估計算法

對于一個平穩(wěn)時間序列y,y(n)表示其在n時刻的值.如果已知時間序列y在某一時刻n之前p個時刻的值，可以建立一個基于該時間序列過去p個值的自回歸模型，并且可用該模型來預測該時間序列將來的值，即

(1)

式中，p為該線性自回歸模型的階數(shù)；ai為y(n-i)項的系數(shù)；e(n)為真實值和預測值之間的誤差項.

若在單變量的基礎上引入另一個平穩(wěn)時間序列x，可建立二維自回歸模型，即

(2)

1.1 自回歸模型的模型階數(shù)估計算法

模型階數(shù)估計是自回歸模型參數(shù)估計過程中的重要步驟.赤池信息準則[7]是常用的模型階數(shù)估計算法，它是一種基于信息論的算法，主要思想是在模型與數(shù)據(jù)的擬合精準度和模型復雜度之間進行權(quán)衡，得出模型的最優(yōu)階數(shù).

給定一個零均值的m維向量自回歸過程zn={zn(1),zn(2),…,zn(m)}T，即

(3)

式中，Φi(i=1,2,…,p)為m×m的系數(shù)矩陣；wn為零均值獨立同分布的向量.假定對于任意的變量i>0，wn和zn-i是相互獨立的.

假設觀測值z1,z2,…,zN都是由式(3)中的自回歸過程所產(chǎn)生，其中N為序列的長度.為了求得模型階數(shù)p的估計值q，設定一個估計階數(shù)的最大上限K，令q取1～K中所有整數(shù)值.對于每一個q值，均采用最小二乘法對原模型進行估計，得到如下模型：

(4)

ZAj=zq+1,jj=1,2,…,m

(5)

式中

式中，Φi(j)表示Φi的第j行.

利用最小二乘法求解式(5)可得

(6)

原模型的預測誤差協(xié)方差矩陣可表示為

(7)

采用AIC定義式選擇最優(yōu)的階數(shù)估計q，即

(8)

令q取遍1～K中所有整數(shù)值，當AIC(q)值最小時，對應的q值即為該模型的最優(yōu)階數(shù).

對于一個二維自回歸模型(見式(2))，由式(8)可知，AIC算法是在假定自回歸模型中不同信號具有相同階數(shù)的基礎上進行估計的，即q=qxx=qyx=qxy=qyy.然而，在實際情況中，自回歸模型中不同信號的階數(shù)往往各不相同，AIC算法會在估計模型階數(shù)時產(chǎn)生過估計，導致后續(xù)系數(shù)估計運算量加大且精準度不高.針對這一問題，本文采用另一種改進的模型階數(shù)估計算法gAIC[13],其定義如下：

(9)

1.2 OPS算法

以式(2)中的時間序列y為例，可得

(10)

采用OPS算法時，首先從式(10)模型項所組成的候選向量中篩選出線性無關項.這些候選向量所構(gòu)成的矩陣表示為

[V1V2…Vqyy+qxy]T

(11)

式(11)中每一行代表一個候選向量.選擇V1作為第1個選出的向量，判定下一個候選向量V2和V1是否線性無關，即利用V1來擬合V2，并計算V2和估計向量的誤差值.一般情況下，在無噪聲的信號中，向量之間應該是線性無關的.但是在有噪聲的情況 (即誤差不為零) 下，可以設置一個最大誤差值，當誤差值小于該最大值時，向量V2便可被選為線性無關的向量，反之則丟棄向量V2.選出V2作為極大線性無關組中的另一個候選向量后，采用向量V1和V2來估計V3的線性無關性.重復以上步驟直到選出一個新的候選向量組合φ=[ω1ω2…ωR]，其中R為所選出的線性無關向量總數(shù).

篩選出新的線性無關的候選向量后，對其進行最小二乘法處理，即將式(10)改寫為

(12)

式中

θg={g1,g2,…,gR}

(13)

式中，gi(i=1,2,…,R)為模型系數(shù).為了最小化誤差ey(t)，令下式達到最小值：

(14)

然后利用最小二乘法可得

(15)

1.3 自適應OPS算法

對于1.2節(jié)中的OPS算法而言，確定候選向量集合，即式(11)中行向量所形成的矩陣，是一個非常重要的過程.AIC準則存在過度擬合的問題，且對二維自回歸模型的每個維度得到相同的階數(shù)估計值，因此，本文采用gAIC算法對自回歸模型進行階數(shù)估計，得到每個維度的各自階數(shù)，繼而進行后續(xù)計算.另一方面，考慮到OPS算法對于噪聲敏感，為了減少噪聲的影響，對原始信號進行重疊加窗處理.以式(10)為例，對于一個長度為N的信號，選擇一個寬度為W的滑動窗，其滑動步長L，具體步驟如下：

①從信號中提取長度為W的信號，記為Wi；

②對每段窗信號Wi，利用gAIC估計出模型階數(shù)，得到式(11)中的矩陣V，進而使用OPS算法計算矩陣V中每個候選項的權(quán)重值；

③將窗向后滑動一段距離L；

⑤計算式(4)中T個窗信號每個候選項權(quán)重的均值，并以此作為該候選項的權(quán)重值，對候選項進行最終篩選.

2 實驗結(jié)果及分析

實驗中采用了如下的二維自回歸模型：

(16)

式中，ex(n)和ey(n)表示均值為0、方差為1的相互獨立的高斯白噪聲序列.

為比較OPS算法和自適應OPS算法的效果，進行了1 000組實驗，并將結(jié)果應用于式(16)y(n)的參數(shù)估計中.在OPS算法中，為了考察模型階數(shù)過大情況下算法的性能，需預設一個比實際情況大的階數(shù)上限.針對本文實驗中所用模型，該階數(shù)上限設為10.由1.2節(jié)可知，模型階數(shù)越大，候選向量越多，從而導致計算量越大.針對這一問題，分別采用AIC和gAIC算法進行階數(shù)估計，以減少不必要的計算.由AIC得到的估計結(jié)果為6，即候選項為Y(n-1),Y(n-2), …,Y(n-6),X(n-1),X(n-2), …,X(n-6)，從而將候選向量個數(shù)由20降至12.采用gAIC算法對模型階數(shù)進行估計，得到的階數(shù)qyy=6，qxy=3，即候選項為Y(n-1),Y(n-2), …,Y(n-6),X(n-1),X(n-2),X(n-3).然后，考慮加窗處理對OPS算法估計模型階數(shù)正確率的影響.由此設計了6種組合實驗方案，見表1.

表1 6種組合實驗方案

下面以表1中第6種組合實驗方案為例，說明模型階數(shù)的選取過程.首先，由gAIC算法得到信號y(n)中的模型階數(shù)qyy=6，qxy=3.然后，由第1節(jié)中所述算法求得每個候選項的權(quán)重值，結(jié)果見圖1.圖中，候選項的權(quán)重值表示每個候選項的重要程度，按遞減順序排列.由圖可知，X(n-3)項的權(quán)重值最大.按照1.2節(jié)所述篩選候選項的原則，權(quán)重值前四大項(即X(n-3)，Y(n-6)，X(n-1)，Y(n-5))相加之和已超出預先設定的閾值0.92，這4項即為最終需保留的候選項.其余5項的權(quán)重值較小，接近于0，表明其對信號Y的貢獻可以忽略不計，該結(jié)果與式(16)中模型一致.

圖1 各候選項的權(quán)重值

為了驗證算法對不同長度信號的有效性，對于表1中的6種方案，每種方案選取2種不同長度(1 000和500)的信號，分別進行1 000次實驗.對于長度為1 000的信號，滑動窗長度為800，步長為40；對于長度為500的信號，滑動窗長度為300，步長為40.考慮到OPS算法需要閾值，而閾值選取沒有統(tǒng)一標準，因此，本文采用不同閾值進行實驗.實驗結(jié)果見圖2.實驗結(jié)果中的正確定義為算法篩選出完全正確的4個候選項，多選和少選都歸為錯誤.由圖可知，大部分閾值范圍(除高閾值范圍)內(nèi)，有窗情況都明顯優(yōu)于無窗情況；在同樣有窗/無窗情況下，gAIC算法的正確率均高于其他2種算法.N=1 000時最佳閾值約為0.92，N=500時最佳閾值約為0.88，正確率均接近90%.當信號長度從1 000減少到500時，OPS算法仍具有較好的性能，如gAIC加窗方案，最佳閾值不同，但正確率并未明顯下降，均接近于90%.

(a) N=1 000

(b) N=500圖2 不同閾值下6種方案結(jié)果對比

圖1中，Y(n-5)項之后候選項的投影值較前4項大幅減少，可以認為候選項集合被分成2個部分，即圖中前4項和后5項.由圖2可知，閾值較大時，算法的正確率急劇降低至0，根據(jù)OPS算法，當閾值慢慢增大直到1時，包含的候選項增加，即逐漸將剩余的權(quán)重值較小的候選項全部選入.由此可知，閾值的選取對算法的正確率起著至關重要的作用.因此，采用了1.3節(jié)中的自適應OPS算法選取最優(yōu)閾值.對于圖1中按權(quán)重值降序排列的9個候選項，依次計算相鄰2項的比值，結(jié)果見圖3.由圖可知，比值為Y(n-5)/X(n-2)時權(quán)重值取得最大值，其后項的權(quán)重值顯著降低.因此，僅保留前面4項.

自適應OPS算法的實驗結(jié)果見圖4.實驗中信號長度為1 000.由圖可知，gAIC有窗的實驗方案取得了最高的正確率，這與圖2(a)中的結(jié)果一致.

圖3 相鄰2項權(quán)重值的比值

綜上可知，通過在算法中引入gAIC，可以有效減少實驗過程中候選項集合里干擾項的個數(shù)；采用加窗的方式，信號能更加平滑，從而減少噪聲帶來的影響.這2種改進均能使OPS算法的正確率明顯提升.自適應OPS算法的正確率較gAIC有窗方案有所下降，但其計算過程中不需要閾值，尤其適用于沒有先驗知識的數(shù)據(jù)或模型.

圖4 自適應OPS算法的實驗結(jié)果

3 結(jié)論

1) 針對自回歸模型參數(shù)估計問題，本文提出了一種基于gAIC和滑動窗的模型參數(shù)估計算法.

2) 在模型參數(shù)估計中，使用gAIC算法對模型階數(shù)進行估計，從而有效減少模型中候選項集合里干擾項的個數(shù).然而，對信號進行加窗處理，使得信號能更加平滑，減少噪聲帶來的影響.最后，采用自適應OPS算法，進一步篩選候選項集合中的候選項，得到最終的模型選項及相應的模型系數(shù)值.

3) 實驗結(jié)果表明，對于2種不同的實驗信號長度，6種組合實驗方案下所提算法都表現(xiàn)出最優(yōu)的魯棒性，其正確率接近于90%.

參考文獻(References)

[1] 蒲慕明, 徐波, 譚鐵牛. 腦科學與類腦研究概述[J]. 中國科學院院刊, 2016, 31(7): 714,725-736. DOI: 10.16418/j.issn.1000-3045.2016.07.001.

[2] Sporns O. Brain connectivity [EB/OL]. (2007-10-27)[2017-05-10]. http://scholarpedia.org/article/Brain_connectivity.

[3] Granger C W J. Investigating causal relations by econometric models and cross-spectral methods [J].Econometrica, 1969,37(3): 424-438. DOI: 10.2307/1912791.

[4] Kamiński M J, Blinowska K J. A new method of the description of the information flow in the brain structures[J].BiologicalCybernetics, 1991,65(3): 203-210. DOI: 10.1007/BF00198091.

[5] Baccalá L A, Sameshima K. Partial directed coherence: A new concept in neural structure determination[J].BiologicalCybernetics, 2001,84(6): 463-474. DOI: 10.1007/PL00007990.

[6] Friston K J. Functional and effective connectivity:A review[J].BrainConnectivity, 2011,1(1): 13-36. DOI: 10.1089/brain.2011.0008.

[7] Akaike H. Power spectrum estimation through autoregressive model fitting[J].AnnalsoftheInstituteofStatisticalMathematics, 1969,21(1): 407-419. DOI: 10.1007/BF02532269.

[8] Rissanen J. A universal prior for integers and estimation by minimum description length[J].TheAnnalsofstatistics, 1983,11(2): 416-431. DOI:10.1214/aos/1176346150.

[9] Khorshidi S, Karimi M, Nematollahi A R. New autoregressive (AR) order selection criteria based on the prediction error estimation[J].SignalProcessing, 2011,91(10): 2359-2370. DOI:10.1016/j.sigpro.2011.04.021.

[10] Zhao Y, Billings S A, Wei H L, et al. A parametric method to measure time-varying linear and nonlinear causality with applications to EEG data[J].IEEETransactionsonBiomedicalEngineering, 2013,60(11): 3141-3148. DOI: 10.1109/TBME.2013.2269766.

[11] Stoica P, Babu P. Model order estimation via penalizing adaptively the likelihood (PAL)[J].SignalProcessing, 2013,93(11): 2865-2871. DOI: 10.1016/j.sigpro.2013.03.014.

[12] Lu S, Ju K H, Chon K H. A new algorithm for linear and nonlinear ARMA model parameter estimation using affine geometry[J].IEEETransactionsonBiomedicalEngineering, 2001,48(10): 1116-1124. DOI: 10.1109/ 10.951514.

[13] Yang C, Le Bouquin Jeannes R, Bellanger J J, et al. A new strategy for model order identification and its application to transfer entropy for EEG signals analysis[J].IEEETransactionsonBiomedicalEngineering, 2013,60(5): 1318-1327. DOI: 10.1109/TBME.2012.2234125.