鄒 岸 林 霖 王 濤

(南方醫(yī)科大學(xué)生物醫(yī)學(xué)工程學(xué)院，廣州，510515)

?

誘發(fā)電位反卷積技術(shù)的不適定問題及正則化解決方法

鄒 岸 林 霖 王 濤

(南方醫(yī)科大學(xué)生物醫(yī)學(xué)工程學(xué)院，廣州，510515)

連續(xù)循環(huán)平均反卷積 (Continuous loop averaging deconvolution, CLAD) 是近年來用于提取高刺激率模式下聽覺誘發(fā)電位(Auditory evoked potential, AEP)的一種行之有效的方法。但是，CLAD方法在頻率域求解時(shí)，對(duì)刺激序列的頻譜特性有嚴(yán)格的限制，給應(yīng)用帶來不便和局限。本文提出一種在時(shí)域?qū)崿F(xiàn)反卷積的方法，將其轉(zhuǎn)化為線性變換矩陣的逆濾波處理。并且利用奇異值分解分析了由不良序列帶來的不適定問題，引入正則化技術(shù)改善病態(tài)矩陣對(duì)重建結(jié)果的影響。最后比較了若干種典型刺激序列和不同噪聲條件下AEP的恢復(fù)實(shí)驗(yàn)，結(jié)果表明本方法可以較好地解決不良序列和一般噪聲水平條件下暫態(tài)AEP信號(hào)的恢復(fù)重建。

聽覺誘發(fā)電位；高刺激率；反卷積；正則化技術(shù)；逆濾波

引 言

當(dāng)聽覺系統(tǒng)受到外界特定的聲音刺激時(shí)，中樞神經(jīng)系統(tǒng)將產(chǎn)生一系列與刺激相關(guān)的微弱生物電反應(yīng)，這種電活動(dòng)可以通過頭皮表面電極記錄將其從腦電背景活動(dòng)中提取出來，稱為聽覺誘發(fā)電位(Auditory evoked potential, AEP)[1]。常規(guī)AEP記錄需要根據(jù)感興趣成分潛伏期的范圍，確定足夠長的刺激間隔(Stimulus onset asynchrony, SOA)，以重復(fù)性短聲作用于聽覺系統(tǒng)，將所記錄的反應(yīng)進(jìn)行疊加平均得到具有暫態(tài)性質(zhì)的AEP。高刺激率AEP(High stimulus rate AEP, HSR-AEP)的記錄條件下，SOA將小于暫態(tài)AEP潛伏期范圍，導(dǎo)致相鄰刺激引發(fā)的反應(yīng)間發(fā)生重迭。這種重迭現(xiàn)象可以看成是暫態(tài)AEP與刺激序列間卷積效應(yīng)的結(jié)果，重建暫態(tài)AEP實(shí)質(zhì)上是對(duì)HSR-AEP進(jìn)行反卷積處理?；谠摾碚摚壳耙延械姆淳矸e算法主要包括最大長序列(Maximum length sequence, MLS)方法[2]、連續(xù)循環(huán)平均反卷積(Continuous loop averaging deconvolution, CLAD)方法[3-4]以及多刺激率穩(wěn)態(tài)平均反卷積(Multi-rate steady-state averaging deconvolution, MSAD)方法[5]。

CLAD采用的是一種非等間隔，即SOA存在抖動(dòng)的刺激序列，將記錄的重迭反應(yīng)在頻域中以逆濾波方式還原暫態(tài)AEP[6-8]。這種做法對(duì)刺激序列的頻譜特性有約束，使得序列的選擇成為影響應(yīng)用的難題[9-11]。本文從CLAD方法最初的時(shí)域模型出發(fā)，在時(shí)域通過線性矩陣逆變換的方式進(jìn)行反卷積重建[3-4]。由于刺激序列構(gòu)造的系統(tǒng)矩陣通常具有不同程度的病態(tài)問題，難以求得一個(gè)穩(wěn)定的最終解,本文提出采用正則化技術(shù)[12]對(duì)求解過程中存在的病態(tài)問題進(jìn)行改善,通過模擬重建實(shí)驗(yàn)對(duì)算法進(jìn)行了驗(yàn)證。

1 CLAD反卷積模型

(1)

2 卷積模型的時(shí)域求解及正則化技術(shù)

H=[hc(t),hc(t-1),hc(t-2),…,hc(t-L)]

(2)

(3)

由式(3)，為了求解AEP只需要在時(shí)域進(jìn)行矩陣逆濾波運(yùn)算

(4)

在實(shí)際應(yīng)用時(shí)，由于采樣率較高，矩陣H將是一個(gè)大型的稀疏矩陣，往往具有不同程度的病態(tài)問題。因此對(duì)于式(4)，在求解過程中需要考慮其不適定性。對(duì)于一個(gè)不適定問題，通常難以求得其精確解。對(duì)于式(4)，將其展開，即

(5)

從式(5)右端第2項(xiàng)可以看出，由于噪聲成分的影響，所求信號(hào)結(jié)果將受到不同程度的擾動(dòng)，在極端情況下，噪聲成分將被嚴(yán)重放大，導(dǎo)致求解結(jié)果失去意義。

針對(duì)上述情況，在對(duì)式(3)求解過程中，本文采用奇異值分解(Singular value decomposition,SVD)及最小二乘思想對(duì)其進(jìn)行分析。對(duì)系統(tǒng)矩陣H進(jìn)行SVD分解如下

H=U∑VT

(6)

(7)

為了減少噪聲影響，求得一個(gè)合理解，可以在對(duì)式(7)的求解中引入正則化技術(shù)，在最小二乘解的基礎(chǔ)上乘以一個(gè)濾波因子fi，減少較小的奇異值對(duì)解的影響。由此，可得正則化解為

(8)

對(duì)式(8)，若忽略奇異值較小項(xiàng)，只考慮前K個(gè)奇異值及其對(duì)應(yīng)奇異向量的作用，即取濾波因子為

(9)

可求得截?cái)郤VD(Truncated SVD,TSVD)解

(10)

(11)

對(duì)應(yīng)正則化解為

(12)

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

3.1CLAD模型仿真

圖1給出本文方法仿真計(jì)算過程和結(jié)果。采用4點(diǎn)B樣條函數(shù)作為基本波形進(jìn)行變換生成理想暫態(tài)中潛伏期AEP信號(hào)，其有效成分由潛伏期及幅值各不相同的5個(gè)主波構(gòu)成，有效長度約為100 ms，模擬采樣率為4 000 Hz(圖1(b))。計(jì)算選用的一組刺激序列總長為248 ms，其刺激所在位置為{0, 32, 62, 98, 132, 158, 186, 224, 248} ms，序列抖動(dòng)率為22.6%，平均刺激率為32.3 Hz(圖1(a))。由該序列構(gòu)成的系統(tǒng)矩陣條件數(shù)為1.4×1018。根據(jù)腦電信號(hào)的性質(zhì)，采用粉紅噪聲生成噪聲信號(hào)，疊加信噪比為0 dB，如圖1(c)所示。暫態(tài)信號(hào)與刺激序列的卷積和噪聲疊加后合成的混合反應(yīng)波形即為仿真所觀測到的信號(hào)如圖1(d)所示。

對(duì)觀測信號(hào)采用不同方式求解的AEP波形如圖1(e)所示。其中細(xì)實(shí)線波形為式(7)的最小二乘解，粗實(shí)線波形為引入正則化后求得的結(jié)果(正則化參數(shù)為0.8)，虛線為真實(shí)解。由于實(shí)驗(yàn)中采用的刺激序列構(gòu)成的系統(tǒng)矩陣具有病態(tài)性，導(dǎo)致直接逆濾波求得的最小二乘解嚴(yán)重偏離真實(shí)值，而引入正則化技術(shù)后對(duì)噪聲擾動(dòng)進(jìn)行了抑制，其解明顯接近于真實(shí)值。

圖1 CLAD模型仿真Fig.1 Simulation for CLAD method

3.2 AEP波形恢復(fù)

圖2 中潛伏期AEP及腦電噪聲Fig.2 Middle-latency AEP and EEG noise

3.2.1 系統(tǒng)矩陣的病態(tài)程度

在CLAD反卷積模型的時(shí)域求解中，系統(tǒng)矩陣的病態(tài)程度直接影響解得質(zhì)量，而該矩陣的構(gòu)成又與刺激序列相關(guān)。本文選用3組具有不同刺激率的刺激序列，由該3組刺激序列構(gòu)成的系統(tǒng)矩陣具有不同的病態(tài)程度，其對(duì)應(yīng)的條件數(shù)和序列的具體參數(shù)見表1,各組總長均為204.8 ms，采樣率為5 000 Hz。

圖3,4,5分別給出了當(dāng)疊加信噪比為0 dB時(shí)，在不同刺激條件下的中潛伏期AEP的恢復(fù)結(jié)果。其中，圖3(a)，4(a)，5(a)表示3組刺激序列，分別對(duì)應(yīng)sq1，sq2及sq3。 圖3(b)，4(b)，5(b)分別表示由sq1，sq2及sq3三組刺激序列構(gòu)成的系統(tǒng)矩陣對(duì)應(yīng)的奇異值取10為底的對(duì)數(shù)后分布情況。在3種刺激條件下，系統(tǒng)矩陣對(duì)應(yīng)都有不同情況的極小奇異值分布，將影響信號(hào)恢復(fù)質(zhì)量。 圖3(c)，4(c)，5(c)分別表示在sq1，sq2，sq3刺激序列條件下，分別采用直接方法求得的最小二乘解和采用正則化技術(shù)求得的結(jié)果對(duì)比，圖中虛線表示原AEP波形，細(xì)實(shí)線表示直接求得的最小二乘解波形，粗實(shí)線表示正則化求得的重建AEP波形。在3種情況下，直接求得的最小二乘解均與真實(shí)解偏差嚴(yán)重，且隨著系統(tǒng)矩陣的條件數(shù)改變其偏差表現(xiàn)也不同。而正則化解雖然與真實(shí)值存在差異，其波形較好地還原了原AEP的基本波形特征，具體表現(xiàn)為在正則化解的波形中，Na，Pa，Nb，Pb等各個(gè)主波均可有效辨識(shí)，恢復(fù)信號(hào)具有較高的可性度。

表1 刺激序列主要參數(shù)

圖3 刺激序列sq1性質(zhì)及其實(shí)驗(yàn)結(jié)果Fig.3 Properties of stimulus sequence sq1 and its experimental results

圖4 刺激序列sq2性質(zhì)及其實(shí)驗(yàn)結(jié)果Fig.4 Properties of sequence sq2 and its experimental results

圖5 刺激序列sq3性質(zhì)及其實(shí)驗(yàn)結(jié)果Fig.5 Properties of stimulus sequence sq3 and its experimental results

3.2.2 噪聲對(duì)重建信號(hào)的影響

針對(duì)刺激序列sq3，分別在3種信噪比情況下采用正則化技術(shù)進(jìn)行反卷積重建暫態(tài)AEP信號(hào)。圖6表示3種噪聲水平下AEP波形的重建結(jié)果，其中圖6(a,c,e)分別為5 dB,0 dB和-5 dB下3種信噪比對(duì)應(yīng)模擬掃程反應(yīng)， 圖6(b，d，f)為上述對(duì)應(yīng)3個(gè)掃程反應(yīng)反卷積后所得的AEP波形。從實(shí)驗(yàn)結(jié)果看，在3種情況下，采用正則化技術(shù)都能在很大程度上抑制噪聲影響，重建所得的AEP在形態(tài)上與原AEP波形接近，各個(gè)主波成分(Na,Pa,Nb,Pb)均能得到辨別。

(13)

從表中數(shù)據(jù)可知，隨著噪聲程度的加重，重建信號(hào)與原信號(hào)間的均方根誤差逐漸變大，同時(shí)相關(guān)系數(shù)及重建信噪比均不同程度的減少。重建信噪比下降幅度較大，提示噪聲對(duì)于恢復(fù)信號(hào)質(zhì)量的影響嚴(yán)重，當(dāng)信噪比減少至-5dB時(shí)，重建信號(hào)雖然在波形特征上依然明顯，但其質(zhì)量要差于前兩者。綜合上述結(jié)果，表明引入正則化技術(shù)后在實(shí)驗(yàn)的3種信噪比下能穩(wěn)定地重建AEP信號(hào)。

圖6 不同信噪比下AEP波形恢復(fù)結(jié)果Fig.6 Deconvolved AEPs under different SNR

表2 不同噪聲下重建AEP信號(hào)與原信號(hào)的對(duì)比

3.3 結(jié)果分析

采用高刺激率方案可以有效減少記錄時(shí)間，也有利于對(duì)聽覺系統(tǒng)進(jìn)行適應(yīng)性研究和病理性診斷[14-15]。CLAD作為高刺激率條件下反卷積重建暫態(tài)AEP信號(hào)的一個(gè)簡便有效的方法，已經(jīng)在許多研究中得到了應(yīng)用[8,16-17]?；跁r(shí)域模型的CLAD反卷積求解，可以轉(zhuǎn)化成線性矩陣的逆濾波處理。根據(jù)文獻(xiàn)[3-4]的理論，采用具有低抖動(dòng)特性的刺激序列可以保證時(shí)域求解中系統(tǒng)矩陣逆的存在性，但其沒有充分考慮噪聲對(duì)恢復(fù)性能的影響。為了避免求逆運(yùn)算過度放大噪聲成分而帶來的失真，需要對(duì)刺激序列的進(jìn)行優(yōu)化選擇。另一方面，良性的逆矩陣往往需要犧牲刺激序列的低抖動(dòng)性要求，或者需要增加一個(gè)掃程的刺激個(gè)數(shù)和刺激長度，這些約束限制了CLAD方法應(yīng)用范圍。本文利用奇異值分解技術(shù)剖析了逆變換矩陣的性質(zhì)，指出刺激序列選擇對(duì)逆矩陣病態(tài)性的影響，并通過引入正則化原理改善病態(tài)矩陣對(duì)信號(hào)重建質(zhì)量的影響，有效地緩解了刺激序列選擇的局限。

本文中AEP波形的恢復(fù)實(shí)驗(yàn)，分別在具有不同病態(tài)程度的系統(tǒng)矩陣和不同噪聲水平下進(jìn)行，其結(jié)果均說明了正則化技術(shù)可以有效抑制矩陣病態(tài)影響。然而，從實(shí)驗(yàn)結(jié)果來看，系統(tǒng)矩陣病態(tài)程度的強(qiáng)弱對(duì)解的影響沒有一個(gè)清晰的判別界限，一個(gè)弱病態(tài)的矩陣在本文實(shí)驗(yàn)中亦有可能導(dǎo)致解的嚴(yán)重失真。這種現(xiàn)象反映了奇異值所影響的成分對(duì)信號(hào)重建的貢獻(xiàn)存在差別，在計(jì)算中必須針對(duì)實(shí)際情況進(jìn)行分析。另外，正則化參數(shù)一方面抑制了噪聲放大，另一方面也影響了信號(hào)的有效成分，其選擇對(duì)求解結(jié)果有著重要影響。采用TSVD和Tikhonov法進(jìn)行正則化求解，其本質(zhì)是統(tǒng)一的。兩者在求解過程中都涉及正則化參數(shù)k，λ的選擇，當(dāng)k值取得偏小或者λ值取得偏大都會(huì)導(dǎo)致過正則化解，反之亦將不能有效抑制噪聲影響。正則化參數(shù)的選取主要包括嶺估計(jì)、L曲線和廣義交叉驗(yàn)證3種方法[18]。然而這些方法確定的正則化參數(shù)并不一致，在實(shí)際應(yīng)用中更有效確定正則化參數(shù)的方法是根據(jù)本文對(duì)待解信號(hào)的先驗(yàn)知識(shí)，比如AEP實(shí)際幅值范圍的了解，結(jié)合上述估計(jì)方法確定正則化參數(shù)的數(shù)值。

4 結(jié)束語

本文研究表明，正則化技術(shù)應(yīng)用于CLAD方法的時(shí)域模型中可以較好地解決信號(hào)重建過程中遇到的病態(tài)問題。運(yùn)用正則化反卷積，可以避免討論序列選擇對(duì)恢復(fù)信號(hào)的影響，當(dāng)系統(tǒng)矩陣具有病態(tài)時(shí)，仍可求得一個(gè)合適解，為在CLAD實(shí)驗(yàn)方式下恢復(fù)暫態(tài)AEP信號(hào)提供了又一個(gè)有效手段。

[1] Kraus N, Nicol T. Auditory evoked potentials[M]. Springer Berlin Heidelberg:Encyclopedia of Neuroscience, 2009:214-218.

[2] Eysholdt U, Schreiner C. Maximum length sequences—A fast method for measuring brainstem-evoked responses [J]. Audiology, 1982, 21(3):242-250.

[3] ?zdamar ?, Delgado R E, Yavuz E, et al. Deconvolution of auditory evoked potentials obtained at high stimulus rates [C]// First International IEEE EMBS Conference on Neural Engineering. Capri Island, Italy： IEEE, 2003:285-288.

[4] Delgado R E, ?zdamar ?. Deconvolution of evoked responses obtained at high stimulus rates [J]. Acoustical Society of America, 2004,115(3):1242-1251.

[5] Wang Tao, Zhan Changan, Yan Gang, et al. A preliminary investigation of the deconvolution of auditory evoked potentials using a session jittering paradigm [J]. Journal of Neural Engineering, 2013, 10(2): 14.

[6] ?zdamar ?, Bohórquez J. Signal-to-noise ratio and frequency analysis of continuous loop averaging deconvolution (CLAD) of overlapping evoked potentials [J]. J Acoust Soc Am, 2009,119(1):429-438.

[7] Wang T, ?zdamar ?, Bohorquez J, et al. Wiener filter deconvolution of overlapping evoked potentials [J]. Neuroscience Methods, 2006, 158(2):260-270.

[8] ?zdamar ?, Bohrquez J, Ray S S. Pb (P1) resonance at 40 Hz: Effects of high stimulus rate on auditory middle latency responses (MLRs) explores using deconvolution [J]. Clinical Neurophysiology, 2007,118(6):1261-1273.

[9] Jewett D L, Caplovitz G, Baird B, et al. The use of QSD (q-sequence deconvolution) to recover superposed, transient evoked-responses [J]. Clinical Neurophysiology, 2004, 115(12): 2754-2775.

[10]Wang T, Huang J H, Lin L, et al. Continuous and discrete-time stimulus sequences for high stimulus rate paradigm in evoked potential studies [J]. Computational and Mathematical Methods in Medicine,2013,2013:1-10.

[11]黃江華, 李彬, 林霖, 等. 改進(jìn)DE算法優(yōu)化高刺激率AEP的刺激序列[J]. 數(shù)據(jù)采集與處理, 2013, 28(5): 672-678.

Huang Jianghua, Li Bin, Lin Lin, et al. Stimulus-sequences optimization for high-rate AEPs using modified differential evolution algorithm [J]. Journal of Data Acquisition and Processing, 2013,28(5):672-678.

[12]Tikhonov A, Arsenin V. Solutions of ill-posed problems [M].Washington DC: Winston, 1977:113-124.

[13]Hansen P C. The discrete picard condition for discrete ill-posed problems [J]. BIT Numerical Mathematics, 1990,30(4):658-672.

[14]Su Yuanyuan, Li Zhenji, Wang Tao. Deconvolution methods and applications of auditory evoked response using high rate stimulation[C]//New Developments in Biomedical Engineering. Croatia:In-Teh, 2010:105-122.

[15]王濤，張愛桃，李彬,等. 希爾伯特-黃變換增強(qiáng)高刺激率AEP的信噪比[J]. 數(shù)據(jù)采集與處理，2011,26(3):280-285.

Wang Tao, Zhang Aitao, Li Bin, et al. Improved signal-to-noise ratio of AEP based on HHT [J]. Journal of Data Acquisition and Processing, 2011,26(3):280-285.

[16]Bohórquez J, ?zdamar ?. Generation of the 40 Hz auditory steady-state response (ASSR) explained using convolution [J]. Clinical Neurophysiology, 2008,119(11):2598-2607.

[17]朱程, 王濤, 黃江華, 等. 40 Hz暫態(tài)聽覺誘發(fā)電位中潛伏期成分的引出分析[J]. 中國生物醫(yī)學(xué)工程學(xué)報(bào), 2013, 32(5):539-545.

Zhu Cheng,Wang Tao, Huang Jianghua, et al. Elieitation analysis for the 40 Hz transient auditory evoked potential [J]. Chinese Journal of Biomedical Engineering, 2013,32(5):539-545.

[18]馮潔婷,顏剛,王濤,等.Tikhonov正則化參數(shù)選擇對(duì)高速率刺激聽覺誘發(fā)電位重建的影響[J]. 航天醫(yī)學(xué)與醫(yī)學(xué)工程,2012,25(1):54-60.

Feng Jieting, Yan Gang, Wang Tao, et al. Effects of parameter selection of Tikhonov regularization on reconstruction of high-rate auditory evoked potentials [J]. Space Medicine & Medical Engineering, 2012, 25(1):54-60.

Regularization Solution to Ill-posedness of Deconvolution Technique for Evoked Potentials

Zou An, Lin Lin, Wang Tao

(School of Biomedical Engineering, Southern Medical University, Guangzhou, 510515, China)

Continuous loop averaging deconvolution (CLAD) is a recently developed method to restore the auditory evoked potential (AEP) under high stimulus rate condition. This method solves the deconvolution problem in frequency domain for computational efficiency, but suffers from stringent limitation in selecting a stimulus sequence with required spectral property. Hereby we propose a new method to solve the deconvolution problem in time domain by constructing a linear transform matrix to model the convolution process. To understand the AEP distortion caused by the ill-posed matrix generated from a bed stimulus sequence, we assess the matrix property using singular value decomposition (SVD) technique and introduce Tikhonov regularization method to deal with the ill-posedness. In the stimulation experiment, we compare some typical sequences with different ill-posedness conditions and restore the transient AEPs under various noise levels. These results justify the proposed approach to the AEP deconvolution with less restriction on the sequence selection.

auditory evoked potential; high stimulus rate; deconvolution; regularization techniques; inverse filtering

國家自然科學(xué)基金(61172033, 61271154)資助項(xiàng)目。

2014-07-24；

2014-09-04

R318

A

鄒岸(1989-)，男，碩士研究生，研究方向：生物電信號(hào)處理、醫(yī)院信息化，E-mail: gdmxza@126.com。

林霖(1981-)，男，講師，研究方向：醫(yī)院信息化、生物電信號(hào)處理。

王濤(1967-)，男，教授，研究方向：生物電信號(hào)處理、神經(jīng)工程。