基于小波分解的顱腦動態電阻抗監測中體動干擾實時處理方法研究

陳曉飛，張 戈，王輝林*

（1.聯勤保障部隊第980 醫院醫學工程科，石家莊 050051；2.聯勤保障部隊第980 醫院放射診斷科，石家莊 050051）

0 引言

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

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

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

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

1 方法

1.1 EIT 的數據采集

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

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

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

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

對于動態EIT 而言，其重構過程可以簡略表示為

式中，y表示阻抗變化圖像；B表示重構矩陣；xf表示前景幀邊界電壓數據；xb表示參考幀邊界電壓數據。體動干擾會體現在前景幀xf中并影響重構結果。

前景幀數據中同時包含了正常的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℃。雙層薄膜覆蓋、露天處理的日較差變化趨勢基本一致，表現為晴天條件下日較差較大，多云次之，陰天最??；而單層薄膜覆蓋表現為多云條件下日較差較大，晴天次之，陰天最小。

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

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

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

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

式中，|wjk|/為標準化后的正態分布；其中x∈R，為正態分布的積分函數。

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

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

圖2 利用小波拒絕概率篩選體動干擾對應的小波分解系數方法示意圖

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

1.3 實時小波計算

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

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

1.4 實驗驗證設計

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

1.4.1 仿真實驗

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

式中，n=4；ω=2πf；μ 為正弦波的振蕩幅度；σ（t）為高斯白噪聲，模擬信號采集過程中的系統噪聲；γ為噪聲幅值，設置為1%水平。xsimulate（t）的幅值控制為（-1，1）。正弦波參數的設置涵蓋可能存在的人體生理信號及其他干擾，具體設置參照文獻[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 個數據點。

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

圖4 仿真信號及體動干擾設置示意圖

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

1.4.2 人體實測數據實驗

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

2 結果

2.1 仿真實驗

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

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

圖5 仿真實驗結果

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

2.2 人體實測數據實驗

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

圖6 人體實測數據處理結果

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

圖7 被試1 電壓均值曲線

圖8 被試2 電壓均值曲線

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

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

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

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

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

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

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

3 討論

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

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

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

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