【作 者】袁思念 ，李若薇 ，朱子孚 ，馬勝才，牛航舵，葉繼倫,，張旭,

1 深圳大學(xué) 醫(yī)學(xué)部 生物醫(yī)學(xué)工程學(xué)院，深圳市，518060

2 深圳市生物醫(yī)學(xué)工程重點(diǎn)實(shí)驗(yàn)室，深圳市，518060

3 廣東省生物醫(yī)學(xué)信息檢測與超聲成像重點(diǎn)實(shí)驗(yàn)室，深圳市，518060

0 引言

全身麻醉是保證患者安全的必要和重要的環(huán)節(jié)，在手術(shù)過程中進(jìn)行準(zhǔn)確的無創(chuàng)麻醉深度監(jiān)測是非常必要的[1]。由于麻醉藥物主要作用于人體的中樞神經(jīng)系統(tǒng)，所以可以通過處理腦電信號(hào)來監(jiān)測患者的麻醉深度，然而腦電信號(hào)十分微弱，并且具有非平穩(wěn)非線性等特性，在實(shí)際手術(shù)中采集到的腦電信號(hào)易受噪聲干擾[2-3]。

在麻醉手術(shù)開始之前，病人此時(shí)為清醒狀態(tài)，容易發(fā)生眨眼等動(dòng)作，引入低頻高幅的眼電信號(hào)，如果不對(duì)眼電信號(hào)進(jìn)行去除，進(jìn)行相應(yīng)的參數(shù)計(jì)算時(shí)低頻段的成分會(huì)明顯偏高；而當(dāng)病人進(jìn)入麻醉狀態(tài)時(shí)，由于手術(shù)原因，采集到的腦電信號(hào)會(huì)有高頻段的肌電以及電刀干擾。由于這些噪聲均是隨機(jī)出現(xiàn)，并且噪聲信號(hào)頻段與腦電信號(hào)的頻段重合，常規(guī)方法不易去除[4]，而麻醉深度監(jiān)測和特征參數(shù)計(jì)算需實(shí)時(shí)進(jìn)行，這些因素給實(shí)時(shí)的噪聲檢測與去除帶來了困難。目前常見的腦電信號(hào)去噪算法有小波閾值法、獨(dú)立成分分析法、自適應(yīng)濾波器等[5]。小波閾值法無法去除眼電以及實(shí)時(shí)檢測噪聲發(fā)生的時(shí)間[6-7]；獨(dú)立成分分析方法要滿足觀測信號(hào)的通道數(shù)大于源信號(hào)數(shù)，并且計(jì)算量較大，會(huì)產(chǎn)生較大的延時(shí)[8-9]；自適應(yīng)濾波器的頻率跟蹤范圍很窄，適用范圍有限[2]。

本研究提出了一種基于離散小波變換自適應(yīng)的提取雙閾值來進(jìn)行噪聲的檢測與去除，該方法能夠?qū)崟r(shí)的檢測出噪聲發(fā)生的時(shí)間并進(jìn)行精確去除，從而提高后續(xù)計(jì)算提取的特征參數(shù)的穩(wěn)定性。

1 方法與步驟

1.1 數(shù)據(jù)預(yù)處理

由于腦電信號(hào)的主要頻帶范圍為0.5～70 Hz，因此要先對(duì)腦電信號(hào)進(jìn)行預(yù)處理。經(jīng)過0.5 Hz高通、70 Hz低通、50 Hz帶阻濾波器，得到原始腦電圖信號(hào)。

1.2 腦電信號(hào)處理流程

圖1是算法的整體流程。將預(yù)處理的腦電信號(hào)經(jīng)過離散小波變換，提取出低頻能量和高頻能量，根據(jù)其所在范圍對(duì)噪聲干擾進(jìn)行分級(jí)，分別對(duì)每種不同的干擾情況進(jìn)行不同的處理，最后將處理后得到的小波系數(shù)進(jìn)行信號(hào)重構(gòu)，所得即為去噪后的信號(hào)，并根據(jù)噪聲干擾的范圍，按照不同的策略對(duì)其閾值自適應(yīng)進(jìn)行實(shí)時(shí)更新。

圖1 算法整體流程Fig.1 The overall flow chart of the algorithm

1.3 離散小波變換

小波變換可以有效地表示信號(hào)的時(shí)頻特征。在低頻時(shí)，小波變換具有較低的時(shí)間分辨率和較高的頻率分辨率，而在高頻時(shí)，小波變換具有較高的時(shí)間分辨率和較低的頻率分辨率[10]。基于Mallat快速算法的離散小波變換是腦電信號(hào)分析的常用算法。利用離散小波變換可以有效地將信號(hào)的頻段與高頻段進(jìn)行分離。根據(jù)下式將腦電信號(hào)分解成細(xì)節(jié)系數(shù)和近似系數(shù)：

這里j是分解級(jí)別，k是對(duì)應(yīng)的小波系數(shù)下標(biāo)；其中，[n]為高通濾波器，將濾掉信號(hào)的低頻分量而輸出高頻分量，h[n]為低通濾波器，濾掉信號(hào)的高頻分量而輸出低頻分量，式中Aj,k與Dj,k分別表示腦電信號(hào)分解得到的近似系數(shù)與細(xì)節(jié)系數(shù)。

從頻率的角度出發(fā)，可以將干擾信號(hào)分為兩類：一類是頻率較低的眼電信號(hào)等造成的干擾，另一類是頻率較高的肌電信號(hào)與電刀工作時(shí)等造成的干擾。

本方法的主要目的是進(jìn)行實(shí)時(shí)噪聲的檢測與去除以提高后續(xù)特征參數(shù)計(jì)算的穩(wěn)定性，為了準(zhǔn)確檢測低頻干擾和高頻干擾發(fā)生的時(shí)間，將數(shù)據(jù)以1 s的長度進(jìn)行分段，1 s的數(shù)據(jù)進(jìn)行離散小波變換，采用db4小波進(jìn)行分解，由于信號(hào)的采樣率為256 Hz，選擇分解層數(shù)為3層，1 s的腦電圖信號(hào)被分解成3組細(xì)節(jié)系數(shù)D1～D3和一組近似系數(shù)A3，其分別對(duì)應(yīng)的頻帶如表1所示。由于腦電信號(hào)主要頻帶范圍為0～70 Hz，所以直接將D1系數(shù)認(rèn)定為噪聲，將其置零。

表1 腦電信號(hào)頻帶分解Tab.1 EEG signal decomposition

1.4 計(jì)算低頻能量與高頻能量

近似系數(shù)A3包含了低頻的能量，細(xì)節(jié)系數(shù)D2和D3則包含高頻能量，通過分解得到的小波系數(shù)，提取這1 s的低頻能量PL與高頻能量PH，計(jì)算公式如下：

這里k是對(duì)應(yīng)的小波系數(shù)下標(biāo)。這樣提取出低頻能量與高頻能量，如果這1 s的信號(hào)受到低頻干擾，則計(jì)算出PL值會(huì)較大，如果有高頻干擾，則計(jì)算出的PH值也會(huì)較大，如果此時(shí)發(fā)生相應(yīng)干擾相當(dāng)于有一個(gè)突變值，因此我們可以檢測出噪聲的發(fā)生。

1.5 噪聲的判斷與去除

為了準(zhǔn)確檢測和去除噪聲，我們使用兩組閾值來判定噪聲，一組閾值用于檢測低頻干擾以及評(píng)定低頻干擾的等級(jí)，定義為TL1與TL2；另一組閾值用于檢測高頻干擾以及評(píng)定高頻干擾的等級(jí)，定義為TH1與TH2。通過使用兩組閾值，提高了噪聲檢測的可靠性，并且減小了因誤判而導(dǎo)致刪除有用信息的概率，閾值的設(shè)定與噪聲的判定規(guī)則如下：

（1）強(qiáng)烈噪聲的判定。

TL2與TH2為強(qiáng)烈噪聲的判定閾值，如果PL＞TL2或PH＞TH2，則判定此時(shí)發(fā)生強(qiáng)烈的干擾，直接去除此段信號(hào)，不參與后續(xù)的參數(shù)計(jì)算。

（2）輕微噪聲的判定。

TL1與TH1為輕微低頻噪聲的判定閾值，如果TL1＜PL＜TL2，則判定此時(shí)發(fā)生輕微的低頻噪聲干擾，此時(shí)，為了盡可能地去除噪聲而保留信號(hào)，如果A3中系數(shù)平方大于閾值TL1的系數(shù)進(jìn)行1/8倍的收縮。同理，如果TH1＜PH＜TH2，則判定此時(shí)發(fā)生輕微的低頻噪聲干擾，此時(shí)，對(duì)D2和D3中系數(shù)平方大于閾值TH1的系數(shù)進(jìn)行1/8倍的收縮。

（3）信號(hào)的判定。

如果PL＜TL1，則判定信號(hào)無低頻干擾，不對(duì)A3的系數(shù)進(jìn)行處理。如果PH＜TH1，此時(shí)判定為信號(hào)沒有高頻干擾，不對(duì)D2與D3的系數(shù)進(jìn)行處理。

經(jīng)過以上的噪聲判定，如果信號(hào)沒有被判定為強(qiáng)烈干擾，則對(duì)經(jīng)過處理后的A3、D2、D3進(jìn)行重構(gòu)，重構(gòu)所得的信號(hào)參與后續(xù)參數(shù)計(jì)算。

1.6 閾值的更新策略

閾值的設(shè)定應(yīng)隨著信號(hào)進(jìn)行自適應(yīng)調(diào)整，閾值TL1與TH1是為了檢測輕微噪聲而設(shè)定的，所以TL1應(yīng)略高于最近一段信號(hào)的平均低頻能量，TH1應(yīng)略高于最近一段信號(hào)的平均高頻能量；而閾值TL2與TH2是為了檢測強(qiáng)烈噪聲而設(shè)定的，TL2應(yīng)較高于最近一段信號(hào)的平均低頻能量，TH2應(yīng)較高于最近一段信號(hào)的平均高頻能量。為了提高閾值計(jì)算的平穩(wěn)性，閾值計(jì)算區(qū)間取30 s。

為了求取最近的30 s信號(hào)的平均低頻能量與平均高頻能量，引入兩組系數(shù)PLtemp與PHtemp，其中，PLtemp為當(dāng)前這1 s中的平均低頻能量，PHtemp為當(dāng)前這1 s中的平均高頻能量；為了對(duì)閾值進(jìn)行平穩(wěn)的更新，根據(jù)當(dāng)前計(jì)算出的PL和PH所在的范圍，設(shè)定了三種更新策略：

（1）策略一。

當(dāng)PL＞TL2時(shí)，此時(shí)判定為強(qiáng)烈低頻干擾，令PLtemp=TLmax；TLmax是能容許的最高的信號(hào)低頻能量，這是一個(gè)經(jīng)驗(yàn)值，一般取200～1000；這里TLmax設(shè)置為500。同理，當(dāng)PH＞TH2時(shí)，此時(shí)判定為強(qiáng)烈高頻干擾，令PHtemp=THmax；THmax是能容許的最高的信號(hào)高頻能量，這也是一個(gè)經(jīng)驗(yàn)值，一般取50～200；這里THmax設(shè)置為100。

（2）策略二。

當(dāng)TL1＜PL＜TL2時(shí)，此時(shí)判定為輕微的低頻干擾，此時(shí)如果PL＞TLmax，令PLtemp=TLmax；否則，PLtemp=PL/2。同理，當(dāng)TH1＜PH＜TH2時(shí)，此時(shí)判定為輕微的高頻干擾，如果PH＞THmax，令PHtemp=THmax；否則，令PHtemp=PH/2。

（3）策略三。

當(dāng)PL＜TL1時(shí)，此時(shí)判定為無低頻干擾，此時(shí)如果PL＞TLmax，令PLtemp=TLmax，否則，令PLtemp=PL。同理，當(dāng)PH＜TH1時(shí)，此時(shí)判定為無高頻干擾，如果PH＞THmax，令PHtemp=THmax；否則，令PHtemp=PH。

數(shù)據(jù)通過不同策略處理后，對(duì)最近30 s的信號(hào)計(jì)算出的30個(gè)PLtemp和30個(gè)PHtemp求平均，分別得到平均低頻能量PLmean和平均高頻能量PHmean；則閾值TL1，TL2，TH1，TH2的計(jì)算分別如式(5)～式(8)所示。

1.7 特征參數(shù)提取

QUAN等[11]認(rèn)為，隨著麻醉深度的增加，腦電信號(hào)的低頻段能量逐漸增加，幅度逐漸增大，高頻能量逐漸減小。BetaRatio是腦電信號(hào)中高頻能量與低頻能量對(duì)數(shù)的比率，能夠衡量腦電信號(hào)能量的變化，GU等[1]認(rèn)為，BetaRatio對(duì)于清醒和麻醉狀態(tài)具有較好的區(qū)分能力。

為了計(jì)算出BetaRatio，首先應(yīng)求其功率譜，為了改進(jìn)譜估計(jì)方差特性，本研究選用的是韋爾奇功率譜估計(jì)方法，窗函數(shù)選擇為漢明窗，計(jì)算長度為10 s，將10 s的數(shù)據(jù)分割成10段，每一段數(shù)據(jù)不重疊，對(duì)每一段數(shù)據(jù)求取功率譜然后求平均，即可得到10 s信號(hào)的功率譜P(k),k=1,2,…,N。

其中，P(k)表示頻率所對(duì)應(yīng)的功率譜值，N為計(jì)算的FFT點(diǎn)數(shù)，這里N=1024；fs為信號(hào)的采樣率，這里fs=256 Hz。

2 結(jié)果

2.1 噪聲的判定與去除結(jié)果

圖2是一段清醒狀態(tài)下的腦電信號(hào)中低頻噪聲的檢測與去除效果，圖3是一段麻醉狀態(tài)下的腦電信號(hào)中高頻噪聲去除效果。從圖2(b)中可以看出，PL在第2秒和第3秒處的值大于TL1但是小于TL2，說明此時(shí)檢測到了這兩秒處發(fā)生了輕微低頻干擾，算法對(duì)這兩秒的低頻系數(shù)進(jìn)行收縮以去除低頻噪聲；圖3(c)中PH在第5秒和第6秒處的值大于TH1但是小于TH2，說明此時(shí)檢測到了這兩秒處發(fā)生了輕微高頻干擾，所以對(duì)這兩秒的高頻系數(shù)進(jìn)行收縮以去除低頻噪聲。圖2(d)和圖3(d)是重構(gòu)后的去噪信號(hào)，從圖中可以看出，噪聲處被精準(zhǔn)去除。

圖2 低頻噪聲的檢測與去除Fig.2 Low frequency noise detection and removal

圖3 高頻噪聲的檢測與去除Fig.3 High frequency noise detection and removal

2.2 參數(shù)的計(jì)算結(jié)果

圖4(a)是使用美國Aspect公司生產(chǎn)的BISX外置模塊及數(shù)據(jù)記錄系統(tǒng)記錄的一例麻醉手術(shù)中采集的腦電信號(hào)，采集過程中BISX模塊濾波器模式設(shè)置為高通0.5 Hz、低通70 Hz、帶阻50和60 Hz，以256 Hz的采樣率記錄。圖4(b)顯示的是計(jì)算出的低頻能量PL與對(duì)應(yīng)的閾值TL1和TL2，圖4(c)顯示的是高頻能量PH與對(duì)應(yīng)的閾值TH1和TH2，從圖中可以看出，兩組閾值根據(jù)信號(hào)的變化自適應(yīng)地調(diào)整大小，檢測出腦電圖信號(hào)中發(fā)生噪聲干擾的位置，并對(duì)受干擾的部分進(jìn)行了相應(yīng)的處理。

圖4 一例病人的腦電信號(hào)去噪處理結(jié)果Fig.4 The denoising result of a patient's EEG signal

BetaRatio與麻醉深度密切相關(guān)[3]，病人較清醒的時(shí)候BetaRatio較高，病人進(jìn)入麻醉狀態(tài)時(shí)，BetaRatio會(huì)進(jìn)入一個(gè)較低的水平，從圖5中可以看出，病人在1000 s左右處進(jìn)入麻醉狀態(tài)，在16000 s左右處開始恢復(fù)意識(shí)。圖5(a)是使用去噪前的數(shù)據(jù)計(jì)算出的BetaRatio，圖5(b)是使用本研究所提出的方法進(jìn)行去噪后所得到的BetaRatio，從圖中可以看出，未經(jīng)去噪的腦電信號(hào)計(jì)算出的BetaRatio在麻醉狀態(tài)下，由于噪聲的干擾，出現(xiàn)很多的異常值，而經(jīng)過本研究提出的自適應(yīng)噪聲檢測與去除算法后，異常值基本被去除，提高了該參數(shù)的穩(wěn)定性。

圖5 參數(shù)BetaRatio去噪前后對(duì)比Fig.5 Comparison of feature parameters BetaRatio before and after denoising

3 結(jié)論

我們提出了一種基于離散小波變換自適應(yīng)的提取雙閾值來進(jìn)行腦電信號(hào)噪聲的檢測與去除方法，將噪聲分為低頻噪聲與高頻噪聲，采用離散小波變換提取對(duì)應(yīng)的低頻能量與高頻能量，并自適應(yīng)地計(jì)算出兩組閾值，采用雙閾值法判定干擾的等級(jí)，提高了噪聲檢測的可靠性；同時(shí)，利用小波變換優(yōu)良的時(shí)頻特性，對(duì)發(fā)生噪聲干擾時(shí)的小波系數(shù)進(jìn)行處理，達(dá)到對(duì)噪聲的精準(zhǔn)去除。最后，對(duì)去噪前后的腦電數(shù)據(jù)進(jìn)行分析計(jì)算特征參數(shù)BetaRatio，結(jié)果表明，去噪后計(jì)算得到的BetaRatio較去噪前異常值明顯減少，參數(shù)的穩(wěn)定性有了進(jìn)一步的提升。