基于高頻中值濾波的小波濾波在腦動脈色素濃度譜特征信號提取中的應(yīng)用*

朱翔宇，姜馨兒，田兆紅，修采，查雨彤

(上海健康醫(yī)學院 上海穿戴式醫(yī)療技術(shù)與器械工程研究中心，上海 201318)

1 引 言

近年來，由于人們生活環(huán)境、生活習慣等因素影響，腦梗塞、腦出血、腦卒中等腦血管疾病和并發(fā)癥的患者明顯增多。在腦血管疾病診斷中，多采用腦血流動力學參數(shù)作為重要的診斷指標。目前，疾病的常規(guī)臨床檢查技術(shù)無論是生化檢測還是影像學檢測，均需要固定在指定區(qū)域，且存在使用成本高，易產(chǎn)生副作用等問題。因此給實時監(jiān)測腦血流動力學參數(shù)增加了難度，耽誤治療。與此同時，診斷的準確性高度依賴于醫(yī)生的診療技術(shù)，增加誤診和漏診的風險。

近紅外光譜技術(shù)(NIRS)與目前常規(guī)的臨床檢測技術(shù)相比，具有瞬時分辨率高、放射性低、無創(chuàng)等優(yōu)點[1]。20世紀80年代末，F(xiàn)errari和Edwards基于菲克原理，提出了利用NIRS測量腦組織能量代謝和血液循環(huán)。這種方法在當時引起了廣大學者的關(guān)注。氧合血紅蛋白(oxyhe myoglobin, HbO2)作為腦血流動力學的示蹤劑，首次用于腦血流動力學的定量分析，但對于腦損傷或其他易受氧飽和度降低影響的疾病患者，腦血流動力學參數(shù)不能對吸收的HbO2的變化做出精確的反饋。

為了解決這些問題，本課題組已經(jīng)設(shè)計了反射式光電傳感測量裝置，并通過提前注射指示劑染料吲哚菁綠(indocyanine green, ICG)能夠進行同步和實時測量[2]。雖然該方法對腦血流動力學參數(shù)的采集十分有效，然而還是包含著組織光吸收干擾、心動、呼吸帶來的干擾，以及一些測量過程中的無關(guān)噪聲。因此，這就需要采取有效的濾波手段提取特征信號[3-4]。目前常見的濾波方法普遍存在濾波后信號易失真，不能有效濾除低頻或者高頻噪聲的問題[5]。為達到理想的濾波效果，本研究采用中值濾波和小波濾波，通過結(jié)合這兩種方法，以實現(xiàn)在同時濾除不同頻率噪聲的基礎(chǔ)上不失真，改善了現(xiàn)有的問題。

2 原理

2.1 腦血流動力學參數(shù)的測量

血流動力學參數(shù)的測量包括腦血流量(CBF),腦血容積(CBV) 和平均通過時間(MTT)等血流參數(shù)。CBF, CBV 和 MTT都是基于菲克原理計算的血流參數(shù)。這些參數(shù)體現(xiàn)了示蹤劑在器官組織中的積累量以及示蹤劑在器官中的到達率與解離率之間的差異[6]。

(1)

其中，Q、CA和CV代表組織、動脈和靜脈中的示蹤劑濃度，F(xiàn)為血流。當t小于通過器官的最小通過時間時，示蹤劑將不會出現(xiàn)在靜脈血液中，故在0到t的時間間隔期內(nèi)靜脈血液中無示蹤劑流出，上式被簡化為式(2)以計算血流量。

(2)

對于注射的示蹤劑，積累量是隨著腦組織中ICG的總濃度ΔC[ICG]tissue增加而變化的，腦內(nèi)ICG動脈引入量表示為ΔC[ICG]blood。根據(jù)菲克定律，CBF以式(3)表達，其中kICG是ICG的分子量。

(3)

CBV是血液與腦組織體積的比值，其計算為測量期間腦組織和動脈血液中ICG濃度積分的比率。

(4)

根據(jù)中心體積原理，平均通過時間(MTT)可以計算如下：

(5)

從上述公式可以看出，腦組織和腦動脈的ICG濃度是腦血流動力學檢測的最關(guān)鍵指標。因此，使用何種方法進行測量，在一定程度上決定了整個腦血流動力學參數(shù)的準確性。通過NIRS技術(shù)，建立腦血流動力學模型，利用光電手段，獲得光電容積脈搏波(PPG)信號，從中可以獲得吸光物質(zhì)的濃度變化情況，進而獲得示蹤劑ICG的濃度。

2.2 NIRS-ICG技術(shù)腦血流動力學模型的建立

腦血管側(cè)支循環(huán)豐富，呈網(wǎng)絡(luò)狀分布，通過檢測動脈的搏動情況能夠獲得大量有價值的信息。脈搏是當大量血液進入動脈后使管徑擴張的節(jié)律性搏動，其光程的變化使血流中光的吸收強度隨之改變。臨床上常利用光電手段獲取PPG用于實時監(jiān)測，測量原理見圖1。

圖1PPG的產(chǎn)生原理圖

Fig.1The generate schematic of PPG

由于頭部組織是一種高散射低吸收的不均勻介質(zhì)，故在進行PPG信號測量時，整個被測區(qū)域可以視為既有光散射又有光吸收能力的對象。因此，吸光物質(zhì)濃度改變所造成的光強變化遵從朗伯-比爾定律[7]。為了修正腦血流動力學狀態(tài)異常條件下CHbT和cSaO2產(chǎn)生的波動，采用776、805和850 nm多波長近紅外光源在測量過程中的時刻，與初始靜息態(tài)的ICG濃度進行差分處理，得到修正后的ICG濃度。

(6)

其中，Ф為光子通量率，T805、T850分別為：

(7)

(8)

3 實驗方法

3.1 小波變換在腦動脈色素濃度譜特征信號濾波中的應(yīng)用

3.1.1小波去噪 對信號進行小波去噪[8]，首先需將原信號進行小波變換。小波變換適合處理短時間內(nèi)發(fā)生突變的信號，能夠充分體現(xiàn)信號的每一個細節(jié)，對于PPG這類信號強度弱，易產(chǎn)生干擾的信號非常適用。利用塔式算法[9]進行小波變換，將一長度為 M 的信號分解為相等長度的小波系數(shù)[10]，并進行后續(xù)的濾波處理。根據(jù)其理論，對腦動脈色素濃度譜特征信號進行小波去噪的步驟以及注意事項如下：

(1)對特征信號做初步分析，采用小波變換[11]對含噪信號進行分解。根據(jù)分解層數(shù)越高，有效信號提取效果越好，均方誤差越小的原則，選擇合適的分解層數(shù)，但考慮到分解層數(shù)過高(大于6)，極易使部分PPG信號失真，致使重構(gòu)后特征信號與原始信號差異過大，影響后續(xù)的參數(shù)分析，所以這里分解層數(shù)應(yīng)等于6；

(2)對于步驟(1)中分解所得到的分解系數(shù)，選擇分層閾值函數(shù)，進行閾值化處理，對各尺度小波系數(shù)去噪。通過分層閾值函數(shù)后的小波系數(shù)能夠在保留信號主要特征的情況下,整體連續(xù)性更佳，便于進行PPG信號重構(gòu)；

(3)對小波進行逆變換，將經(jīng)分層閾值函數(shù)處理過的小波系數(shù)加以重構(gòu)，由此得到小波濾波去噪后的特征信號。

3.1.2小波濾波中閾值函數(shù)的選擇 在閾值函數(shù)的選擇中，有兩個較為常用的閾值函數(shù)，軟閾值和硬閾值，但均不能滿足該信號去噪的要求，對于該含噪信號，在進行去噪處理后，既要使去噪后信號的信噪比較高，還要考慮所得信號的平滑性。因此，在此引入了分層閾值函數(shù)，該函數(shù)結(jié)合了軟閾值和硬閾值函數(shù)去噪的優(yōu)點，在保留了盡可能多的信號特征的情況下，使重構(gòu)后的信號過渡更加平滑，其表達式為：

(9)

式中，σ為噪聲標準方差，j為小波分解層數(shù)，Tj為閾值。這里的N根據(jù)上文有關(guān)分解層數(shù)值的分析，取6為宜。這種方式減少了閾值處理的固定偏差，更好地保留了脈搏波信號在分解各層的特征，有效提升了降噪效果。

3.2 中值濾波

中值濾波是一種基于排序統(tǒng)計理論的能有效抑制噪聲，并對脈沖信號有良好濾除作用的非線性信號處理技術(shù)。這種濾波方式的基本原理是將數(shù)字圖像或數(shù)字序列中一點的值，用該點的一個鄰域中各點值的中值來代替[12]，使周圍的像素值接近真實值，由此達到消除孤立噪聲點的效果。

中值濾波的數(shù)學表示為：

(10)

3.3 結(jié)合中值濾波和小波濾波的高頻中值小波濾波

小波濾波在濾波處理上有獨特優(yōu)勢，能對成分進行獨立分析，但不能對粗差的剔除做到完善。中值濾波在有效解決粗差的剔除上，具有獨特優(yōu)勢[13]。因此，本研究采用了基于高頻中值濾波的小波濾波，用于處理采集的PPG信號，改善了僅采用小波濾波的局限性。

首先在改進算法中，對PPG信號用中值濾波處理，得到平滑且已經(jīng)剔除高頻噪聲的信號。然后將處理后的信號進行小波分解，得到相對應(yīng)的分解系數(shù)，再通過恰當?shù)拈撝岛瘮?shù)[14]進行處理。最后將得到的系數(shù)重構(gòu)，得到所需信號。

根據(jù)上述理論，可將基于高頻中值濾波的小波濾波流程進行歸納，見圖2。

圖2高頻中值濾波的小波濾波步驟

Fig.2The steps of wavelet filtering of high frequency median filter

具體步驟：

(1)對含有噪聲的信號作初步分析，由于信號中的噪聲主要以高頻系數(shù)形式出現(xiàn)，所以，使用中值濾波處理信號中的高頻系數(shù)部分；

(2)將得到的高頻系數(shù)與信號中的低頻系數(shù)加以重構(gòu)，得到一個完成初步濾波的信號；

(3)在步驟(2)得到的信號上再使用小波濾波進行處理；

(4)選取合適的小波基，確定小波分解系數(shù)；

(5)對于步驟(4)中分解所得的分解系數(shù)，選擇合適的軟閾值或硬閾值進行閾值化處理，對各尺度小波系數(shù)去噪[15]；

(6)對小波進行逆變換，將經(jīng)軟硬閾值處理過的小波系數(shù)加以重構(gòu)，最終得到基于中值濾波的小波濾波的特征信號。

4 實驗過程

4.1 實驗的評價手段

為了驗證所提出方法的有效性，本研究利用表征觀察值與真值離散均方根誤差。其中，RMSE越小，信號處理結(jié)果越接近目標信號。

(11)

式中，N是樣本xn的個數(shù)。

同時利用所得的均方根誤差，計算信噪比，信噪比越大，說明混在信號里的噪聲越小。信噪比單位為分貝(dB)。

(12)

式中，Ps和Pn分別代表信號和噪聲的有效功率。

信號的頻譜圖表明了信號在不同頻率分量成分的大小，相比時域圖像提供了更豐富具體的頻域圖像?？梢岳秒x散傅立葉變換實現(xiàn)時快速算法[16]即快速傅里葉變換(fast fourier transform, FFT)，對所得信號進行頻譜分析，在要求精度不高的情況下，可在頻譜圖中讀出信號的成分，觀察頻譜圖可以看出高頻信號是否在經(jīng)過中值濾波后被有效剔除。

(13)

式中，F(xiàn)(ω)為f(t)的像函數(shù)，f(t)為F(ω)的像原函數(shù)。

4.2 仿真實驗

為了更好地比較不同濾波方式所達到的濾波效果，在使用從人體獲取的真實脈搏波進行數(shù)據(jù)處理之前，利用功能強大的MATLAB，使用仿真技術(shù)構(gòu)建調(diào)試環(huán)境，對真實信號進行仿真[17-18]。通過仿真技術(shù)模擬實驗所獲得的參數(shù)結(jié)果，為真實數(shù)據(jù)的處理提供必要參考。

4.2.1信號仿真的步驟 本研究按照PPG信號的大致趨勢以及真實信號噪聲成分復雜的特征進行構(gòu)建，使構(gòu)建完畢的仿真信號與真實信號相比能夠獲得理想的均方根誤差值以及相關(guān)系數(shù)。

為了還原采樣時實測光電容積脈搏波信號在一個周期內(nèi)的變化，在構(gòu)建時，選擇了將正弦信號與余弦信號相疊加[19]，以還原真實信號走勢，見圖3(a)。

由于PPG信號是心臟跳動所引發(fā)的波動信號，因此，該信號的主要頻率首先應(yīng)根據(jù)心臟搏動所產(chǎn)生的心率來決定。臨床上，我們往往用脈率(pulse rate，PR)來替代，根據(jù)成人在靜息狀態(tài)PR參考值為60次/min，構(gòu)建脈搏容積波頻率為1 Hz的仿真信號，見圖3(b)。

由于實測信號在測量中受到多種干擾因素影響，尤其是呼吸產(chǎn)生的干擾[20]，進而使用高斯噪聲來還原這種干擾。高斯噪聲是一種較為泛用的噪聲模型，根據(jù)成人的呼吸頻率以及呼吸噪聲[21]的不規(guī)則性分布，在信號不同的區(qū)間上，分別疊加了均值為0，標準差分別為30、20、5、15、19、2的正態(tài)分布高斯噪聲，見圖3(c)。

疊加后，圖像呈現(xiàn)了較為明顯的不規(guī)則性。此外在測量過程中還存在工頻干擾以及測量設(shè)備傳輸過程中所產(chǎn)生的干擾，根據(jù)其特征，疊加了伽馬噪聲，該伽馬噪聲由三個服從指數(shù)分布的噪聲疊加而成，對此仿真圖像的構(gòu)建具有十分顯著的效果,見圖3(d)。

仿真完成后，對構(gòu)建完畢的仿真信號進行評價，將仿真信號與通過實測獲得的真實信號進行對比，經(jīng)計算，得到RMSE為2.38，相關(guān)系數(shù)為0.928，符合實驗要求。

圖3PPG仿真信號疊加

Fig.3Superposition of simulation PPG signals

4.2.2特征信號的提取 在MATLAB中導入了構(gòu)建的仿真信號。為將其干擾噪聲有效剔除，選擇了常見且實用的中值濾波。將仿真信號放入中值濾波器。中值濾波后的圖像見圖4。

圖4中值濾波效果圖

Fig.4Median filter rendering

從圖中可以看出，中值濾波對高頻噪聲有著較為顯著的濾波效果。對比原信號和濾波后的頻譜圖，在高頻區(qū)域內(nèi)的噪聲已經(jīng)被完全剔除，符合實驗目的。但是根據(jù)中值濾波后呈現(xiàn)出的頻譜圖像，見圖5，能夠發(fā)現(xiàn)低頻信號并未被有效過濾，如呼吸、心動等頻率低的干擾信號仍對特征信號的提取有一定影響。因此在中值濾波后再次選用了小波濾波。

圖5 原信號與中值濾波后的信號頻譜圖

利用小波濾波進行實驗時，首先需要選用合適的小波基和分解系數(shù)(N=6)，之后對中值濾波后的信號進行分解，可以得到6層分解信號和位于第一層的中值濾波后信號，分解信號的圖像按照頻率由低到高依次排列。在閾值作用方式上，采用了能使圖像更平滑的軟閾值函數(shù)這一作用方式，同時也避免了硬閾值“一刀切”帶來的影響。在閾值的選取方面，使用了分層閾值函數(shù)進行閾值選取操作，能使濾波后信號的相似程度更高。根據(jù)式(11)計算可得，每一層閾值分別為99.6031、74.6282、62.7461、55.4213、50.3208、46.5034。這樣能使所得信號在相似度較高的基礎(chǔ)上，確保光滑度。

通過觀察中值濾波和組合濾波后信號的頻譜分析圖，可以看出，之前未被有效剔除的低頻干擾信號，在結(jié)合小波濾波后被有效去除，效果見圖6。綜上，此組合方法能有效濾除仿真光電容積脈搏波的干擾噪聲，可嘗試將其運用于真實的光電容積脈搏波信號濾波。

圖6中值濾波后與組合濾波后的信號頻譜圖對比(低頻部分)

Fig.6Comparison of signal spectrum after median filtering and combined filtering (low frequency part)

4.3 實測信號實驗

按照仿真實驗所采取的過程和方法，本研究利用仿真信號所得的參數(shù)用于實測信號。實驗結(jié)果及頻譜分析見圖7。采用五種方法對特征信號提取后的RMSE、SNR和r見表1。通過分析圖、表可知實驗結(jié)果較為理想，說明該特征信號的噪聲頻段分布廣，采用傳統(tǒng)的低通濾波及高通濾波難以有效濾除。另外，采用傳統(tǒng)的小波濾波方法，選取通用閾值難以有效濾除不同頻段的噪聲所固有的一些特征，易造成偏差，導致失真。而基于高頻中值濾波的小波濾波能夠有效濾除各種因素造成的干擾噪聲，在重構(gòu)后，原始信號的局部特征被很好地保留下來，失真更小。此組合濾波方式對于腦動脈色素濃度譜特征信號的提取是有效的。

表1 種方法特征信號提取后的RMSE、SNR和r

圖7組合濾波效果圖及頻譜分析

Fig.7Combined filter effect and spectrum analysis

5 總結(jié)

基于NIRS-ICG色素濃度譜技術(shù)用于腦血流動力學參數(shù)測量，雖能建立光電容積脈搏波的模型，但其使用的常規(guī)濾波方法提取特征信號仍存在高頻信號去噪能力弱、低頻信號無法準確濾除的問題。為了避免上述問題，本研究提出的基于高頻中值濾波的小波濾波方式，將兩種已經(jīng)比較完善的除噪方法結(jié)合使用，從最終的結(jié)果看，效果優(yōu)于使用單個方法，具有輸出信噪比高、均方差低的優(yōu)勢。從頻譜分析圖來看，各個頻段的干擾都被有效濾除。但是，本研究提取出的特征信號，在低頻部分仍有一部分失真有待改善。在后續(xù)的研究中，要找出低頻信號失真的原因，并加以改進。本研究提供了一種新的濾波思路，為下一步的研究提供方向，希望能給其他學者提供些許參考。