基于小波包分解的腦部EIT測量干擾分析與抑制方法研究

梁 文，胡世頡，劉學(xué)超，徐燦華，付 峰，楊 濱*

（1.空軍軍醫(yī)大學(xué)軍事生物醫(yī)學(xué)工程學(xué)系，西安 710032；2.西安電子科技大學(xué)生命科學(xué)技術(shù)學(xué)院，西安 710126；3.空軍軍醫(yī)大學(xué)西京醫(yī)院神經(jīng)外科，西安 710032）

0 引言

電阻抗斷層成像（electrical impedance tomography，EIT）是一種以體內(nèi)組織生理或病理性改變引起的阻抗變化為目標(biāo)的醫(yī)學(xué)成像技術(shù)，其具有無輻射損傷、時間分辨力高、設(shè)備性價比高、可進行連續(xù)實時成像的優(yōu)勢[1]，目前已應(yīng)用于肺部[2]、腦部[3]、腹部[4]等器官功能和損傷的圖像監(jiān)測中。然而，在腦部EIT的臨床監(jiān)護中，頭部激勵-測量電極經(jīng)常受到體動干擾的影響，導(dǎo)致測量信號失真。體動干擾主要來源于醫(yī)護人員的治療、護理操作和患者有意識或者無意識的體動，這都會引起電極接觸阻抗的變化，在EIT中表現(xiàn)為明顯的偽影，嚴重降低了圖像表征病情變化的能力[5]。

針對上述顱腦EIT的體動干擾，既往研究基于加權(quán)相關(guān)系數(shù)法對單一脫落電極進行檢測，并且使用灰色模型預(yù)測方法，對電極脫落后的數(shù)據(jù)進行補償，以恢復(fù)應(yīng)有的監(jiān)護圖像；同時，使用主成分分析（principal component analysis，PCA）方法和因子矯正法對不同情況下的體動干擾信號進行預(yù)處理[6]。在此基礎(chǔ)上，最近的一項研究使用循環(huán)策略實現(xiàn)多個脫落電極的實時檢測，并初步分析了體動干擾對解調(diào)后一維數(shù)據(jù)的影響，使用曲線擬合差分的方法來抑制體動干擾造成的基線變化，將受體動影響期間的數(shù)據(jù)進行還原處理并優(yōu)化了數(shù)據(jù)質(zhì)量[7-8]。

根據(jù)EIT的測量原理，用于成像的邊界電壓數(shù)據(jù)是由采集到的激勵響應(yīng)信號經(jīng)數(shù)字解調(diào)而來。然而，上述方法對于臨床顱腦EIT體動干擾的處理均是針對解調(diào)后的電壓數(shù)據(jù)，只分析測量信號保留的幅值信息，而忽略了原始信號的頻率特性以及在解調(diào)過程中存在的干擾。本研究將基于顱腦EIT原始測量電壓信號的解調(diào)過程，分析原始測量信號中體動干擾的特征，設(shè)計消除體動干擾信號的預(yù)處理算法，并對算法進行仿真以及臨床試驗驗證。

1 EIT原始電壓信號的測量方法

根據(jù)腦組織阻抗頻譜特性，顱腦EIT常使用1~200 kHz正弦恒流信號作為激勵[9]，該信號注入頭部后在頭皮表面形成相應(yīng)正弦電壓信號。在16電極系統(tǒng)下，選擇一對電極作為激勵，同時測量其余各相鄰電極對上的電壓信號。當(dāng)激勵電極對遍歷所有電極后，可得到192組電壓信號（對向電極對激勵、鄰近電極對測量）。EIT系統(tǒng)的測量模塊使用模數(shù)轉(zhuǎn)換器（analog-to-digital converter，ADC）對每組電壓信號進行每周期128個點的采樣，采樣數(shù)據(jù)經(jīng)正交序列數(shù)字解調(diào)得到用于顱腦EIT圖像重構(gòu)的幅值和相位數(shù)據(jù)。當(dāng)系統(tǒng)工作頻率較高時，則采用基于異步欠采樣的等效時間采樣策略對電壓信號進行離散化[10]，以降低對高分辨力ADC的性能需求。該策略的原理為：設(shè)定采樣信號在一定周期內(nèi)保持不變，信號周期為Ts，2個相鄰采樣點的間隔時間為Tt，每周期需要完成N個點的采樣，現(xiàn)以Tt=13×Ts/N的時間間隔進行采樣。當(dāng)N=128時，則經(jīng)過13個信號周期后可完成128個不相重復(fù)的點的采樣，而且第129個點又回到第1個點的相位。這些點與在同一信號周期內(nèi)以Tt'=Ts/N的間隔采樣結(jié)果相一致，因而可以對信號進行重建，等效時間采樣策略示意圖如圖1所示。由此可以使用26個信號周期完成2個周期256個點的采樣任務(wù)。但是體動干擾也會對等效時間采樣過程造成影響，因此既往研究僅對解調(diào)后的測量數(shù)據(jù)進行處理是片面的。

圖1 等效時間采樣策略示意圖

2 顱腦EIT原始測量信號干擾分析和處理方法

2.1 顱腦EIT系統(tǒng)和測量數(shù)據(jù)

本研究的臨床數(shù)據(jù)來自課題組自研的腦部高精度EIT數(shù)據(jù)采集系統(tǒng)[9]，該系統(tǒng)可工作在1~190 kHz的頻率區(qū)間，輸出電流在0~1.25 mA間程控可調(diào)，使用16電極系統(tǒng)進行顱腦組織阻抗的測量。原始邊界電壓信號通過一個16 bit采樣精度和1 MS/s采樣速率的ADC進行采集，數(shù)據(jù)采集策略為基于異步欠采樣的等效時間采樣。在臨床數(shù)據(jù)的采集中，將電極帶配置于被測者眉弓和枕骨粗隆以上1~2 cm的環(huán)面上并用彈力繃帶固定（如圖2所示），激勵電流頻率為50 kHz，1 s內(nèi)完成每幀196組邊界電壓信號的采集。EIT監(jiān)測開始后，分別記錄因護理操作（如翻身、抬頭等）造成被測志愿者體動干擾的數(shù)據(jù)，待后續(xù)分析和處理。

圖2 電極系統(tǒng)配置效果圖

2.2 體動干擾數(shù)據(jù)的判別方法

根據(jù)原始電壓信號的特點，將標(biāo)準正弦信號與原始電壓信號進行皮爾遜相關(guān)分析，將相關(guān)系數(shù)作為判斷測量信號是否為干擾信號的參考，其中皮爾遜相關(guān)系數(shù)是衡量隨機變量X與Y相關(guān)程度的一種方法，相關(guān)系數(shù)的取值范圍為[-1，1]，其絕對值越大，則表明X與Y相關(guān)度越高[11]。當(dāng)X與Y呈線性相關(guān)時，相關(guān)系數(shù)ρXY取值為1（正線性相關(guān)）或-1（負線性相關(guān)），其計算公式如下：

式中，Cov（X，Y）為X、Y的協(xié)方差；E（·）為括號內(nèi)數(shù)據(jù)的數(shù)學(xué)期望；D（·）為括號內(nèi)數(shù)據(jù)的方差為括號內(nèi)數(shù)據(jù)的標(biāo)準差。

顱腦EIT算法是根據(jù)前景幀與背景幀的差分進行重構(gòu)成像，因此根據(jù)背景幀中邊界電壓的解調(diào)幅值A(chǔ)0j，每幀信號的數(shù)據(jù)集S的計算公式如下：

式中，k為EIT采集的幀數(shù)；j為一幀中192個一維曲線值；ω為激勵電流的頻率；i為每個通道數(shù)據(jù)的采樣值（i=1，2，…，256）。

根據(jù)公式（1），前景幀原始電壓信號與各標(biāo)準信號的相關(guān)系數(shù)集Nj以及解調(diào)值差分百分比集Dj的計算公式如下：

式中，j為前景幀的個數(shù)，j=1，2，3；Sj為第j個前景幀的原始數(shù)據(jù)集；S0為模擬的標(biāo)準信號數(shù)據(jù)集；Aj為第j個前景幀解調(diào)值的集；A0為背景幀解調(diào)值的集。根據(jù)Nj、Dj的分析確定區(qū)分正常邊界電壓信號和體動干擾信號的相關(guān)系數(shù)閾值fth。

2.3 受干擾原始電壓信號的特征分析和小波包分解

分別對正常和受干擾的原始電壓信號進行傅里葉變換，分析主頻信號的頻譜特征，識別干擾信號所處的頻段。利用小波包分解[10，12]剔除干擾頻段，重構(gòu)正常頻段的信號，達到濾除干擾信號的目的。其中小波包分解又稱為最優(yōu)子帶樹結(jié)構(gòu)，是對小波變換的進一步優(yōu)化[13]，其結(jié)構(gòu)如圖3所示。小波包分解的算法思想是在小波變換的基礎(chǔ)上，在每一級信號分解時，除了對低頻子段進行進一步分解，也對高頻子段進行進一步分解。然后通過一個代價函數(shù)計算出最優(yōu)的信號分解路徑，并以此分解路徑對原始信號進行分解。

圖3 小波包分解結(jié)構(gòu)

可將尺度函數(shù)φ（t）記為v0（t），小波函數(shù)ψ（t）改記為v1（t），并設(shè)濾波器的系數(shù)分別為h（k）、g（k），其中h（k）滿足∑h（k-2n）h（k-2m）=δn，m（δn，m為脈沖信號，n與m根據(jù)所選的小波包種類確認），∑h（k）=對于原來的二尺度關(guān)系，可以重新定義遞歸函數(shù)如下：

其中，v0（t）與v1（t）滿足二尺度伸縮方程：v1（t）=∑kgkv0（2t-k）。

小波包分解采用Mallat算法，根據(jù)公式（5）~（6），采用如下遞歸方程實現(xiàn)：

式中，W0為遞歸方程的初始項；與分別為分解層數(shù)為j的偶數(shù)項與奇數(shù)項；x（i）為原始信號，i=0，1，2，…，N/2-1，n=0，1，2，…，2j-1-1，其中N為采樣點數(shù)，j為分解層數(shù)。

根據(jù)上述的遞歸方程可以將信號正交分解到相對應(yīng)的頻段上。小波包分解算法如下：

根據(jù)小波包分解與重構(gòu)原理，設(shè)計受干擾信號的識別與預(yù)處理算法流程如下：

（1）受干擾信號的識別。計算原始電壓信號與標(biāo)準正弦信號之間的相關(guān)系數(shù)a，根據(jù)fth判斷是否對信號進行處理。

（2）基于小波包分解的第一步進行預(yù)處理。采用db8小波對原始電壓信號進行5層分解，共分為32（25）個頻段。每個頻段的帶寬為16 Hz，由于標(biāo)準的原始電壓信號頻率在0~16 Hz范圍內(nèi)，使用節(jié)點為（0，16）、帶寬范圍為0~16 Hz的小波包進行重構(gòu)。針對處理之后的結(jié)果，截取中間2個周期的信號作為最終第一步預(yù)處理結(jié)果。

（3）對第一步預(yù)處理結(jié)果進行分析、反饋。計算經(jīng)過第一步預(yù)處理后的信號與標(biāo)準正弦信號的相關(guān)系數(shù)b，根據(jù)fth判斷是否需要對第一步處理的結(jié)果進行二次處理。

（4）針對反饋結(jié)果進行最終干擾抑制處理。當(dāng)反饋參數(shù)b的值小于fth，則判斷該組數(shù)據(jù)失效，并使用背景幀中的數(shù)據(jù)對該數(shù)據(jù)進行補充。

3 試驗結(jié)果及分析

3.1 臨床干擾數(shù)據(jù)判別閾值的確定

本研究使用課題組自研的型號為FMMU-EIT5的腦部高精度EIT數(shù)據(jù)采集系統(tǒng)實時測量EIT數(shù)據(jù)。選取5名健康志愿者分別編號為1~5，使其呈仰臥位進行約5 min的顱腦EIT動態(tài)監(jiān)測。其中1~4號志愿者在監(jiān)測過程中無體動行為，分別采集630幀原始電壓信號，5號志愿者首先采集50幀無體動行為時的原始電壓信號，然后使其進行翻身和抬頭動作，模擬臨床護理操作過程，再采集286幀原始電壓信號。試驗共采集5組共計2 856幀原始電壓信號，其中第1~4組以及第5組中前50幀數(shù)據(jù)為正常測量數(shù)據(jù)，第5組后286幀數(shù)據(jù)為受體動干擾測量數(shù)據(jù)。

根據(jù)成像結(jié)果，由受體動干擾的原始測量信號解調(diào)得到的數(shù)據(jù)會造成EIT圖像有明顯的偽影。如圖4所示，在相同窗寬下圖4（a）因體動干擾出現(xiàn)嚴重的失真，而圖4（b）無明顯的偽影。

圖4 正常與干擾情況下的重構(gòu)圖像結(jié)果

在5號志愿者的數(shù)據(jù)中分別選擇3幀正常原始測量數(shù)據(jù)和1幀受體動干擾的原始測量數(shù)據(jù)進行分析，計算相關(guān)系數(shù)集Nj以及解調(diào)值差分百分比集Dj，將Dj中所有數(shù)據(jù)作為x軸、Nj中所有數(shù)據(jù)作為y軸，畫出散點圖，如圖5所示。其中，圖5（a）為整體結(jié)果圖，5（b）為x軸在[-10，10]區(qū)間的局部結(jié)果圖。由所有志愿者數(shù)據(jù)的分析結(jié)果得出絕大部分的解調(diào)值差分百分比集中在[-10，10]，而解調(diào)值差分百分比與相關(guān)系數(shù)之間有一定的關(guān)系。因此，將Dj分為集合[-10，10]和[-∞，-10）∪（10，∞]，并歸納出正常原始測量數(shù)據(jù)和受體動干擾的原始測量數(shù)據(jù)的相關(guān)系數(shù)分布，詳見表1。同時求取5組數(shù)據(jù)中正常原始電壓信號相關(guān)系數(shù)的平均值與方差，詳見表2。

圖5 相關(guān)系數(shù)與解調(diào)值差分百分比散點圖

根據(jù)表1的結(jié)果以及正常測量數(shù)據(jù)幀與受干擾數(shù)據(jù)幀成像結(jié)果可以分析得出：當(dāng)相關(guān)系數(shù)達到0.996時，重構(gòu)圖像就可以較好地消除偽影。同時根據(jù)表2的結(jié)果可以分析得出，正常原始電壓信號的相關(guān)系數(shù)集中分布在0.996周圍。因此，將相關(guān)系數(shù)0.996定義為區(qū)分正常信號和體動干擾信號的閾值fth。

表1 相關(guān)系數(shù)的分布統(tǒng)計單位：%

表2 5組數(shù)據(jù)中正常原始電壓信號相關(guān)系數(shù)的統(tǒng)計學(xué)結(jié)果

3.2 體動干擾的抑制

（1）體動干擾信號的頻譜分析。

對受干擾數(shù)據(jù)幀中相關(guān)系數(shù)低于閾值fth的原始電壓信號進行頻譜分析。結(jié)果表明，與標(biāo)準激勵正弦信號相比，等效時間采樣的電壓數(shù)據(jù)在波形和頻譜上均存在明顯干擾。正常與受體動干擾的原始電壓信號時、頻域特征如圖6所示。因此，利用小波包分解方法來區(qū)分有用頻段與干擾頻段，并對低頻段的信號進行算法重構(gòu)，可達到消除干擾的目的。

圖6 正常與受體動干擾的原始電壓信號時、頻域結(jié)果

（2）小波包分解消除體動干擾。

使用小波包分解方法對志愿者試驗中判定為受體動干擾的原始電壓信號進行處理。從存在干擾的數(shù)據(jù)幀中任意選擇4組（分為A、B、C、D組）受干擾原始電壓信號，對比其處理前、后的時域結(jié)果圖如圖7所示，其相關(guān)系數(shù)以及解調(diào)值差分百分比結(jié)果詳見表3。從表3中可以看出，小波包分解和重建方法使得原始電壓信號與標(biāo)準正弦信號的相關(guān)系數(shù)得到提升，更接近預(yù)設(shè)的閾值fth。同時解調(diào)值差分百分比的絕對值較處理前大幅度降低，即受干擾原始電壓信號的幅值更接近于前景幀原始電壓信號的幅值，而這有效地抑制了等效時間采樣電壓信號中的體動干擾。

表3 4組受干擾原始電壓信號處理前后的相關(guān)系數(shù)以及解調(diào)值差分百分比結(jié)果

圖7 4組受干擾原始電壓信號處理前后的時域結(jié)果

將所有相關(guān)系數(shù)低于閾值fth的采集數(shù)據(jù)進行小波包分解處理后，原始電壓信號的相關(guān)系數(shù)分布得到明顯改善，詳見表4。對于正常數(shù)據(jù)來說，相關(guān)系數(shù)大于fth（0.996）的占比由60.68%提高到72.00%，0.950~0.996的占比由14.32%降低到3.00%。受體動干擾數(shù)據(jù)的相關(guān)系數(shù)則全部≥0.950，其中＞0.996的占比由18.75%提高到24.61%，僅有0.39%處于0.950~0.996范圍。本研究中的小波包分解方法能有效改善顱腦EIT原始電壓測量數(shù)據(jù)的質(zhì)量，尤其是能抑制體動干擾對測量數(shù)據(jù)的影響。

表4 處理后原始電壓信號相關(guān)系數(shù)的分布統(tǒng)計單位：%

將處理之后的數(shù)據(jù)進行正交序列解調(diào)處理，并通過EIT阻尼最小二乘成像算法進行圖像重構(gòu)，結(jié)果如圖8所示。對比處理前的成像結(jié)果（如圖4所示），可以明顯地看出體動干擾在EIT圖像中產(chǎn)生的偽影得到了有效抑制。

圖8 處理后的重構(gòu)圖像

4 結(jié)語

本研究針對顱腦EIT臨床使用中存在的采集數(shù)據(jù)受體動干擾問題，通過分析正常與受體動干擾的原始電壓測量數(shù)據(jù)的差異，提出使用皮爾遜相關(guān)系數(shù)判斷受干擾數(shù)據(jù)的方法，并基于小波包分解與重建方法實現(xiàn)對體動所產(chǎn)生測量數(shù)據(jù)偏置性干擾的抑制作用。志愿者監(jiān)測試驗結(jié)果表明，本研究所提方法能夠有效地修正EIT原始電壓信號等效時間采樣中存在的體動干擾，改善采樣信號的頻譜質(zhì)量，提高EIT效果。由于臨床護理操作和患者體動會使EIT圖像產(chǎn)生嚴重偽影，影響了顱腦EIT對患者病情變化的實時監(jiān)測，因此消除和抑制體動干擾、提升信號采樣精準度對顱腦EIT的臨床應(yīng)用具有重要意義。此外，EIT電極和皮膚接觸阻抗的變化是導(dǎo)致測量數(shù)據(jù)體動干擾的根本原因，而本研究所提方法并不能評估EIT電極與皮膚接觸阻抗的變化量，因此下一步將結(jié)合電極-皮膚接觸阻抗的測量方法，進一步優(yōu)化EIT原始測量信號干擾的識別與抑制方法。