基于小波分解的顱腦動態(tài)電阻抗監(jiān)測中體動干擾實時處理方法研究

陳曉飛，張 戈，王輝林*

（1.聯(lián)勤保障部隊第980 醫(yī)院醫(yī)學(xué)工程科，石家莊 050051；2.聯(lián)勤保障部隊第980 醫(yī)院放射診斷科，石家莊 050051）

0 引言

顱腦動態(tài)電阻抗斷層成像（electrical impedance tomography，EIT）是一種無創(chuàng)、低成本、可連續(xù)監(jiān)護的生物醫(yī)學(xué)功能成像技術(shù)，在腦損傷等疾病的早期診斷方面有著良好的應(yīng)用前景[1-2]。動態(tài)EIT 系統(tǒng)通過均勻安放在頭部的16 個電極，在2 個不同的時刻向顱腦施加安全電流激勵并測量邊界電壓，并利用一定的算法根據(jù)邊界電壓重構(gòu)出2 個時刻間顱內(nèi)的阻抗變化分布[3]。為了保證重構(gòu)圖像的真實有效，參考幀時刻和前景幀時刻的邊界電壓差應(yīng)盡可能只來自于顱內(nèi)阻抗變化，而良好穩(wěn)定的數(shù)據(jù)采集過程是滿足這一需求的首要條件[4]。然而，在長時連續(xù)監(jiān)護的過程中，電極易受臨床上的體動因素如患者的體動（如擺頭、翻身）以及醫(yī)護人員的護理操作等的影響，使得電極-皮膚的接觸狀態(tài)發(fā)生改變，從而在邊界電壓的測量中引入體動干擾，影響重構(gòu)圖像的質(zhì)量[5-7]。因此，需要對信號中的體動干擾進行處理。

體動干擾是人體生理相關(guān)光、電信號采集過程中的常見問題。以是否需要額外的參考信號輸入為標準，可將體動干擾處理方法大致分為2 類：第一類為基于自適應(yīng)濾波器的處理方法[8-9]。其特點在于需要額外的輸入信號作為信號處理的期望，以此作為體動干擾或理想信號的先驗信息再對干擾信號進行處理。此類方法要求硬件設(shè)備具備相應(yīng)的參考信號采集功能，因此其適用范圍受到一定限制。第二類處理方法主要包括維納濾波法[10]、信號相關(guān)提升法[11]、卡爾曼濾波法[12]、主成分分析法[13]、獨立成分分析法[14]等。以上方法不需要額外的輸入信號，而是通過利用信號的時域、空域或頻域特征，或采用奇異值分解等方法，對體動干擾信號進行檢測、分離。然而，上述方法均有其特定的使用前提，并且難以進行實時處理[15]。

小波分解法是一種通過分析信號的時頻特性來檢測并處理體動干擾成分的信號處理方法。最早Sato 等[16]利用小波分解的方法實現(xiàn)了顱腦近紅外光光譜檢測中體動偽影的精準識別。后續(xù)研究人員在此基礎(chǔ)上通過對對應(yīng)體動干擾的小波分解系數(shù)進行處理，實現(xiàn)體動干擾的抑制[9，17-18]。Yang 等[19]通過小波分解的方法對肺部呼吸EIT 監(jiān)測中體動干擾信號進行處理，表明小波分解法在阻抗監(jiān)測信號的質(zhì)量提升方面也有著良好的應(yīng)用潛力。目前，將小波分解方法應(yīng)用到顱腦電阻抗監(jiān)測中的文章鮮有報道。

基于此，本研究利用小波分解的方法處理顱腦EIT 監(jiān)測中的體動干擾信號。與傳統(tǒng)小波信號處理的方法不同，本研究擬基于不同來源信號小波分解系數(shù)的分布特征對體動干擾信號進行檢測，并采用壘墻式的計算策略實現(xiàn)該處理方法的實時運行。本方法擺脫了傳統(tǒng)小波處理方法在軟硬閾值選擇上的局限性，通過概率分布的方法實現(xiàn)小波分解系數(shù)的動態(tài)選擇，提高了數(shù)據(jù)處理的效能。

1 方法

1.1 EIT 的數(shù)據(jù)采集

本研究使用的顱腦EIT 系統(tǒng)需要等間距受試者頭部環(huán)貼16 個電極，采用對向激勵、鄰近采集的工作模式進行邊界電壓數(shù)據(jù)的測量，數(shù)據(jù)采集流程如圖1 所示。該數(shù)據(jù)采集模式下單次激勵可得到12 個邊界電壓數(shù)據(jù)，遍歷所有電極構(gòu)成完整的一幀測量數(shù)據(jù)[20]。

圖1 對向激勵-鄰近測量模式下EIT 數(shù)據(jù)采集示意圖

1.2 基于小波分解的體動干擾處理

通常顱腦EIT 監(jiān)測采集得到多個通道平緩慢變的時間序列邊界電壓信號[21]。當(dāng)體動干擾發(fā)生時，邊界電壓信號發(fā)生劇烈波動，與腦電等人體生理相關(guān)電信號中體動干擾信號的形態(tài)相似，使得使用小波分解的方法來區(qū)分正常EIT 信號與體動干擾成為可能。

對于動態(tài)EIT 而言，其重構(gòu)過程可以簡略表示為

式中，y表示阻抗變化圖像；B表示重構(gòu)矩陣；xf表示前景幀邊界電壓數(shù)據(jù)；xb表示參考幀邊界電壓數(shù)據(jù)。體動干擾會體現(xiàn)在前景幀xf中并影響重構(gòu)結(jié)果。

前景幀數(shù)據(jù)中同時包含了正常的EIT 信號和體動干擾信號，可以表示為

式中，xf（t）表示時間序列的前景幀邊界電壓信號；xf0（t）表示正常的EIT 邊界電壓信號；xf0（t）表示正常的EIT 信號；ε（t）表示體動干擾信號。使用Mallet 方法進行小波分解的形式可以將原始信號表示為[22-23]

從圖3可以看出，試驗期間，不同覆蓋條件下，水溫平均日較差從高到低為單層薄膜覆蓋、對照、雙層薄膜覆蓋；單層薄膜覆蓋條件下，30、60 cm水深的池子水溫平均日較差分別為5.9、5.8℃， 相差不大；雙層薄膜覆蓋條件下，30、60 cm水深的池子水溫平均日較差分別為2.5、2.2℃，同樣相差不大；但露天條件下，60 cm 水池的變化范圍在 2.3～7.0℃，平均日較差為4.3℃，比30 cm的水池高1.1℃。雙層薄膜覆蓋、露天處理的日較差變化趨勢基本一致，表現(xiàn)為晴天條件下日較差較大，多云次之，陰天最??；而單層薄膜覆蓋表現(xiàn)為多云條件下日較差較大，晴天次之，陰天最小。

式中，?jk（t）=2j/2?（2jt-k），表示重構(gòu)尺度函數(shù)；ψjk（t）=2j/2ψ（2jt-k），表示重構(gòu)小波函數(shù)；j和k分別表示小波分解層數(shù)和尺度平移系數(shù)；v表示尺度分解系數(shù)；wj表示第j層小波分解系數(shù)，并且有

式中，g（n）和h（n）分別為高通濾波器和低通濾波器。

任意時間序列信號的小波分解系數(shù)均可以使用包含2 個獨立的高斯分布模型的混合高斯分布模型來描述其分布特征[24-25]。2 個高斯分布模型的均值均為0 但方差不同。方差較小的高斯模型代表幅值變化幅度較小的小波分解系數(shù)分布特征，方差較大的高斯模型則代表幅值變化幅度較大的小波分解系數(shù)分布特征。與體動干擾信號相比，正常的顱腦EIT 信號在時間序列幅值改變更為緩慢，因此其對應(yīng)的小波分解系數(shù)方差更小，在0 均值周圍波動范圍更窄。所以，可使用小方差的單個高斯分布模型來描述正常顱腦EIT 信號小波分解系數(shù)的分布[24]。對應(yīng)wjk的分布可以表示為wjk~N（0，σ2j）。其中每層的方差σ2j估計依據(jù)絕對中位差確定。絕對中位差對數(shù)據(jù)集中的異常值不敏感，具有良好的穩(wěn)定性。在顱腦EIT 數(shù)據(jù)采集過程中，體動干擾表現(xiàn)為疊加在正常顱腦EIT 信號上的離散異常值，其對應(yīng)的總體分布方差與期望估計的先驗正態(tài)分布方差不同。少量體動干擾的發(fā)生不影響對正常顱腦EIT 信號小波分解系數(shù)先驗高斯分布σ2j的估計結(jié)果。σ2j的具體計算公式為[26]

每個wjk對應(yīng)一個曲線下面積pjk，將其稱為拒絕概率，表示wjk并非來源于預(yù)定義正態(tài)分布總體的概率，具體定義為

式中，|wjk|/為標準化后的正態(tài)分布；其中x∈R，為正態(tài)分布的積分函數(shù)。

參考統(tǒng)計置信區(qū)間定義，設(shè)定置信概率閾值α，若pjk≥α，則表明wjk并非來源于預(yù)定義正態(tài)分布總體，其對應(yīng)的原始數(shù)據(jù)為體動干擾信號，否則對應(yīng)正常的顱腦EIT 信號，基于此對wjk進行處理，具體可表示為

包含體動干擾的wjk來源于先驗正態(tài)分布N的事件發(fā)生概率位于置信區(qū)間外，其對應(yīng)的拒絕概率應(yīng)大于置信概率閾值，因此受體動干擾影響嚴重的分解系數(shù)置零，從而抑制體動干擾成分。閾值α決定了對體動干擾的判定標準和抑制程度，當(dāng)α→0 時，則對所有的小波分解系數(shù)均不作處理。本研究中設(shè)置α=0.9。各參數(shù)定義如圖2 所示。

圖2 利用小波拒絕概率篩選體動干擾對應(yīng)的小波分解系數(shù)方法示意圖

對小波分解系數(shù)的處理與否取決于體動干擾的嚴重程度。體動干擾的嚴重程度具體定義為當(dāng)前尺度下小波分解系數(shù)超出預(yù)定義閾值α的數(shù)量。定義數(shù)據(jù)集合和，其中I為指示函數(shù)。設(shè)置尺度系數(shù)j，得到數(shù)據(jù)集合{ψj}，進行體動干擾篩選處理的小波分解層標記為ψj，processed，{ψj，processed}的和不應(yīng)小于{ψj}的90%，以保證體動干擾處理的效果，并據(jù)此反向確定小波分解的層數(shù)。

1.3 實時小波計算

通常離散小波分解采用的是Mallet 方法，其通過濾波器族以逐級遞推的方式進行計算。在進行離線小波分解計算時，一般將全部數(shù)據(jù)一次性讀入后再逐層分解。然而，由于分解和重構(gòu)濾波器的長度有限，單次卷積計算所需要的數(shù)據(jù)長度也是有限的，數(shù)據(jù)長度滿足向下一層分解時即可進行小波分解運算?；诖?，饒貴安等[27]提出了一種壘墻式的小波分解計算方法（如圖3 所示），隨著數(shù)據(jù)的不斷讀入，同時進行小波分解的計算，數(shù)據(jù)長度滿足計算需求即向下進行一次小波分解計算，計算流程整體橫向發(fā)展。

圖3 壘墻式小波分解示意圖

1.4 實驗驗證設(shè)計

為了驗證上述方法的有效性，開展仿真實驗和人體實測數(shù)據(jù)實驗研究。

1.4.1 仿真實驗

設(shè)計以下含高斯白噪的復(fù)合頻率正弦信號模擬顱腦EIT 信號：

式中，n=4；ω=2πf；μ 為正弦波的振蕩幅度；σ（t）為高斯白噪聲，模擬信號采集過程中的系統(tǒng)噪聲；γ為噪聲幅值，設(shè)置為1%水平。xsimulate（t）的幅值控制為（-1，1）。正弦波參數(shù)的設(shè)置涵蓋可能存在的人體生理信號及其他干擾，具體設(shè)置參照文獻[28]，其中包括：（1）同步心率信號，f1=1 Hz，μ1=0.4；（2）同步呼吸信號，f2=0.25 Hz，μ2=0.4；（3）低頻混疊信號，f3=0.1 Hz，μ3=0.4；（4）極低頻混疊信號，f4=0.01 Hz，μ4=0.5。仿真信號共包含1 000 個數(shù)據(jù)點。

向上述仿真信號中添加尖峰信號模擬體動干擾。本研究設(shè)置如圖4 所示的4 種不同類型的尖峰信號，分別為MA1、MA2、MA3、MA4。4 種尖峰信號與原始信號疊加得到混合信號xmix。采用均方根誤差百分比（percent root difference，PRD）、Pearson 相關(guān)系數(shù)（r）、決定系數(shù)（R2）3 個量化指標分別對混合信號x和處理后信號y與原始信號xsimulate的一致性進行評價，量化評估模擬體動干擾的處理效果。其中PRD 反映2 組數(shù)據(jù)的離差特征，PRD 越小則2 組數(shù)據(jù)的背離程度越??；r評價2 組數(shù)據(jù)的相似程度，r越大則2 組數(shù)據(jù)越相似；R2反映在空間距離上y（t）向x（t）的逼近程度，最大值為1，越接近1 說明y（t）對x（t）的擬合程度越好。具體的計算公式如下：

圖4 仿真信號及體動干擾設(shè)置示意圖

式中，x（t）和y（t）為要比較的時間序列信號，N為數(shù)據(jù)長度。對于本研究，x（t）為不含模擬體動干擾的原始信號xsimulate（t），y（t）為處理前后的含體動干擾信號xmix（t）。

1.4.2 人體實測數(shù)據(jù)實驗

被試的納入標準如下：（1）頭皮表面無創(chuàng)傷；（2）非危重癥，無顱骨骨折；（3）意識清醒，可配合進行數(shù)據(jù)采集；（4）無其他不適合進行顱腦EIT 數(shù)據(jù)采集的情況。最終納入的被試為2018 年3 月于空軍軍醫(yī)大學(xué)附屬醫(yī)院神經(jīng)外科就診的2 例男性患者（年齡分別為18 歲、53 歲）的數(shù)據(jù)，均為輕型頭部外傷后留觀。將16 個進行過消毒處理的杯狀電極涂抹導(dǎo)電膏后，均勻環(huán)貼于患者頭部。數(shù)據(jù)采集使用FMMU-EIT5 數(shù)據(jù)采集系統(tǒng)，激勵電流為500 μA，激勵頻率為50 kHz，數(shù)據(jù)采集速度為1 幀/s。在數(shù)據(jù)采集開始前，將16 個進行過消毒處理的杯狀電極涂抹導(dǎo)電膏后，均勻環(huán)貼在患者頭部。數(shù)據(jù)采集過程中患者躺在病床上，對患者的體動和醫(yī)護人員的護理操作不做干涉。人體數(shù)據(jù)采集已獲倫理委員會批準。

2 結(jié)果

2.1 仿真實驗

本研究中數(shù)據(jù)處理方法對仿真信號的處理效果如圖5 所示。仿真實驗采用db6 小波，分解層數(shù)為10。圖5（a）、（b）、（c）、（d）為圖4 中對應(yīng)的混合信號經(jīng)過小波處理后的結(jié)果。在信號形態(tài)上，尖峰信號模擬的體動干擾得到了有效抑制，處理后的數(shù)據(jù)與原始數(shù)據(jù)高度相似。量化評估體動干擾抑制方法的結(jié)果見表1。不同指標的改善程度與干擾信號類型密切相關(guān)。對于MA1和MA2，ΔR2大幅上升，而對于MA3和MA4，則Δr改善最明顯。除MA2以外，ΔPRD 的改善均超過70%。

表1 混合信號處理前后PRD、r 及R2 量化指標的比較

圖5 仿真實驗結(jié)果

表1 為混合信號處理前后均方根誤差百分比PRD、Pearson 相關(guān)系數(shù)r、決定系數(shù)R2量化指標的比較。ΔPRD、Δr、ΔR2為對應(yīng)指標在數(shù)據(jù)處理后相較于處理前的百分比改變。

2.2 人體實測數(shù)據(jù)實驗

人體實測數(shù)據(jù)的處理結(jié)果如圖6 所示。由于顱腦EIT 有192 個通道的有效數(shù)據(jù)，因此需要對192個通道分別進行處理。針對納入的2 名被試，分別截取650 幀和280 幀長度的臨床實測數(shù)據(jù)，其中包含了數(shù)個觀測到的體動干擾。選用db4 小波，分解層數(shù)為6，對數(shù)據(jù)進行處理，處理結(jié)果如圖6 所示。對比處理前后的信號波形，原始數(shù)據(jù)中包含的體動干擾都得到抑制，數(shù)據(jù)的連續(xù)性明顯提升。

圖6 人體實測數(shù)據(jù)處理結(jié)果

顱腦EIT 圖像監(jiān)護需要將邊界電壓數(shù)據(jù)重構(gòu)為二維阻抗分布圖像，因此對處理前后的數(shù)據(jù)進行圖像重構(gòu)以檢驗干擾抑制效果。圖像重構(gòu)使用基于Matlab 2016b 的EIDORS v3.8 工具箱，圖像重構(gòu)方法為高斯牛頓法。圖7、圖8 分別為被試1、被試2 的電壓均值曲線圖。圖9~12 為被試1 數(shù)據(jù)處理前后的圖像重構(gòu)結(jié)果對比，圖13~15 為被試2 數(shù)據(jù)處理前后的圖像重構(gòu)結(jié)果對比。體動干擾導(dǎo)致重構(gòu)圖像出現(xiàn)明顯的體動偽影[如圖9（a）、10（a）、11（a）、12（a）所示]，經(jīng)小波方法處理后圖像的均勻性顯著提升[如圖9（b）、10（b）、11（b）、12（b）所示]。在相同坐標尺度范圍下，數(shù)據(jù)處理后的重構(gòu)圖像三維幅值分布更加平滑。同樣的，從重構(gòu)圖像重構(gòu)值變化范圍上看，在相同尺度下，處理后數(shù)據(jù)的重構(gòu)圖像上的偽影基本消失，重構(gòu)圖像的幅值更為平滑。

圖7 被試1 電壓均值曲線

圖8 被試2 電壓均值曲線

圖9 被試1 體動干擾M1 處理前后的重構(gòu)圖像及其三維幅值分布

圖10 被試1 體動干擾M2 處理前后的重構(gòu)圖像及其三維幅值分布

圖11 被試1 體動干擾M3 處理前后的重構(gòu)圖像及其三維幅值分布

圖12 被試1 體動干擾M4 處理前后的重構(gòu)圖像及其三維幅值分布

圖13 被試2 體動干擾M1 處理前后的重構(gòu)圖像及其三維幅值分布

圖14 被試2 體動干擾M2 處理前后的重構(gòu)圖像及其三維幅值分布

圖15 被試2 體動干擾M3 處理前后的重構(gòu)圖像及其三維幅值分布

3 討論

本研究采用小波分解的方法對顱腦EIT 數(shù)據(jù)采集中的尖峰型體動干擾進行處理。與既往小波處理方法不同的是，本研究放棄了既往研究中需要手動確定硬閾值或軟閾值的方法，采用了概率計算的策略進行體動干擾信號的篩選?；谡ｏB腦EIT 信號的小波分解系數(shù)構(gòu)建代表其總體的高斯分布模型，針對單個小波分解系數(shù)，計算其在先驗高斯模型下對應(yīng)的曲線下面積，即為拒絕概率。如果該小波分解系數(shù)位于預(yù)設(shè)定的先驗高斯分布置信區(qū)間外，則其對應(yīng)的拒絕概率應(yīng)大于預(yù)設(shè)的概率閾值，由此實現(xiàn)體動干擾成分的篩選并加以處理。該方法僅需設(shè)定統(tǒng)一的概率閾值，不需要為每層小波分解系數(shù)單獨設(shè)置特定閾值，大大減輕了算法設(shè)計的負擔(dān)，提高了體動干擾的檢測效能。特別是對于EIT 信號這樣的多通道數(shù)據(jù)而言，該方法尤其適用。人體實測數(shù)據(jù)實驗表明，該方法可有效抑制時間序列信號中的體動干擾成分，恢復(fù)信號的連續(xù)性，減輕重構(gòu)圖像偽影。

前期文獻對于體動干擾信號的形態(tài)特征進行過經(jīng)驗性總結(jié)，研究人員將生理相關(guān)光電信號采集中的體動干擾總結(jié)為時間序列信號中幅值超出正常信號范圍的快速、劇烈變化信號成分，并且與正常信號相比更缺乏周期性[16]。在顱腦EIT 信號中，與平緩慢變的正常信號相比，體動干擾主要包含高頻成分。因此，根據(jù)小波函數(shù)與信號間斷點的卷積特性，體動干擾部分信號對應(yīng)的小波分解系數(shù)表現(xiàn)為局部的模極大值[19]。這種小波分解系數(shù)模值的差異是使用混合高斯分布模型描述小波分解系數(shù)分布特征并對干擾信號進行鑒別的基礎(chǔ)。

本研究仍存在一定的不足：本研究中提出的基于小波分解的尖峰型體動干擾處理方法在干擾界定的閾值上具有較好的靈活性，但其大范圍推廣應(yīng)用仍需進行更多的臨床實測數(shù)據(jù)驗證。此外，臨床上體動干擾的發(fā)生有可能存在尖峰型干擾與其他類型體動干擾混疊的場景，降低本研究中小波方法的處理效能。因此，在下一步的研究中將收集更多的臨床體動干擾類型數(shù)據(jù)，對本研究中處理方法進行針對性的改進，增強魯棒性，從而提升顱腦EIT 的臨床實用性與適用性。

綜上，本研究提出的基于小波分解系數(shù)統(tǒng)計特征的顱腦EIT 體動干擾處理方法避免了對每次小波分解結(jié)果設(shè)定具體閾值，而是通過設(shè)定置信概率區(qū)間的方法實現(xiàn)體動干擾成分的篩選與處理，特別適合顱腦EIT 多通道數(shù)據(jù)采集的應(yīng)用場景。實驗結(jié)果表明，該方法可有效抑制體動干擾成分，恢復(fù)正常的圖像監(jiān)護。