朱翔宇,姜馨兒,田兆紅,修采,查雨彤
(上海健康醫學院 上海穿戴式醫療技術與器械工程研究中心,上海 201318)
近年來,由于人們生活環境、生活習慣等因素影響,腦梗塞、腦出血、腦卒中等腦血管疾病和并發癥的患者明顯增多。在腦血管疾病診斷中,多采用腦血流動力學參數作為重要的診斷指標。目前,疾病的常規臨床檢查技術無論是生化檢測還是影像學檢測,均需要固定在指定區域,且存在使用成本高,易產生副作用等問題。因此給實時監測腦血流動力學參數增加了難度,耽誤治療。與此同時,診斷的準確性高度依賴于醫生的診療技術,增加誤診和漏診的風險。
近紅外光譜技術(NIRS)與目前常規的臨床檢測技術相比,具有瞬時分辨率高、放射性低、無創等優點[1]。20世紀80年代末,Ferrari和Edwards基于菲克原理,提出了利用NIRS測量腦組織能量代謝和血液循環。這種方法在當時引起了廣大學者的關注。氧合血紅蛋白(oxyhe myoglobin, HbO2)作為腦血流動力學的示蹤劑,首次用于腦血流動力學的定量分析,但對于腦損傷或其他易受氧飽和度降低影響的疾病患者,腦血流動力學參數不能對吸收的HbO2的變化做出精確的反饋。
為了解決這些問題,本課題組已經設計了反射式光電傳感測量裝置,并通過提前注射指示劑染料吲哚菁綠(indocyanine green, ICG)能夠進行同步和實時測量[2]。雖然該方法對腦血流動力學參數的采集十分有效,然而還是包含著組織光吸收干擾、心動、呼吸帶來的干擾,以及一些測量過程中的無關噪聲。因此,這就需要采取有效的濾波手段提取特征信號[3-4]。目前常見的濾波方法普遍存在濾波后信號易失真,不能有效濾除低頻或者高頻噪聲的問題[5]。為達到理想的濾波效果,本研究采用中值濾波和小波濾波,通過結合這兩種方法,以實現在同時濾除不同頻率噪聲的基礎上不失真,改善了現有的問題。
血流動力學參數的測量包括腦血流量(CBF),腦血容積(CBV) 和平均通過時間(MTT)等血流參數。CBF, CBV 和 MTT都是基于菲克原理計算的血流參數。這些參數體現了示蹤劑在器官組織中的積累量以及示蹤劑在器官中的到達率與解離率之間的差異[6]。
(1)
其中,Q、CA和CV代表組織、動脈和靜脈中的示蹤劑濃度,F為血流。當t小于通過器官的最小通過時間時,示蹤劑將不會出現在靜脈血液中,故在0到t的時間間隔期內靜脈血液中無示蹤劑流出,上式被簡化為式(2)以計算血流量。
(2)
對于注射的示蹤劑,積累量是隨著腦組織中ICG的總濃度ΔC[ICG]tissue增加而變化的,腦內ICG動脈引入量表示為ΔC[ICG]blood。根據菲克定律,CBF以式(3)表達,其中kICG是ICG的分子量。
(3)
CBV是血液與腦組織體積的比值,其計算為測量期間腦組織和動脈血液中ICG濃度積分的比率。
(4)
根據中心體積原理,平均通過時間(MTT)可以計算如下:
(5)
從上述公式可以看出,腦組織和腦動脈的ICG濃度是腦血流動力學檢測的最關鍵指標。因此,使用何種方法進行測量,在一定程度上決定了整個腦血流動力學參數的準確性。通過NIRS技術,建立腦血流動力學模型,利用光電手段,獲得光電容積脈搏波(PPG)信號,從中可以獲得吸光物質的濃度變化情況,進而獲得示蹤劑ICG的濃度。
腦血管側支循環豐富,呈網絡狀分布,通過檢測動脈的搏動情況能夠獲得大量有價值的信息。脈搏是當大量血液進入動脈后使管徑擴張的節律性搏動,其光程的變化使血流中光的吸收強度隨之改變。臨床上常利用光電手段獲取PPG用于實時監測,測量原理見圖1。

圖1PPG的產生原理圖
Fig.1The generate schematic of PPG
由于頭部組織是一種高散射低吸收的不均勻介質,故在進行PPG信號測量時,整個被測區域可以視為既有光散射又有光吸收能力的對象。因此,吸光物質濃度改變所造成的光強變化遵從朗伯-比爾定律[7]。為了修正腦血流動力學狀態異常條件下CHbT和cSaO2產生的波動,采用776、805和850 nm多波長近紅外光源在測量過程中的時刻,與初始靜息態的ICG濃度進行差分處理,得到修正后的ICG濃度。
(6)
其中,Ф為光子通量率,T805、T850分別為:
(7)
(8)
3.1.1小波去噪 對信號進行小波去噪[8],首先需將原信號進行小波變換。小波變換適合處理短時間內發生突變的信號,能夠充分體現信號的每一個細節,對于PPG這類信號強度弱,易產生干擾的信號非常適用。利用塔式算法[9]進行小波變換,將一長度為 M 的信號分解為相等長度的小波系數[10],并進行后續的濾波處理。根據其理論,對腦動脈色素濃度譜特征信號進行小波去噪的步驟以及注意事項如下:
(1)對特征信號做初步分析,采用小波變換[11]對含噪信號進行分解。根據分解層數越高,有效信號提取效果越好,均方誤差越小的原則,選擇合適的分解層數,但考慮到分解層數過高(大于6),極易使部分PPG信號失真,致使重構后特征信號與原始信號差異過大,影響后續的參數分析,所以這里分解層數應等于6;
(2)對于步驟(1)中分解所得到的分解系數,選擇分層閾值函數,進行閾值化處理,對各尺度小波系數去噪。通過分層閾值函數后的小波系數能夠在保留信號主要特征的情況下,整體連續性更佳,便于進行PPG信號重構;
(3)對小波進行逆變換,將經分層閾值函數處理過的小波系數加以重構,由此得到小波濾波去噪后的特征信號。
3.1.2小波濾波中閾值函數的選擇 在閾值函數的選擇中,有兩個較為常用的閾值函數,軟閾值和硬閾值,但均不能滿足該信號去噪的要求,對于該含噪信號,在進行去噪處理后,既要使去噪后信號的信噪比較高,還要考慮所得信號的平滑性。因此,在此引入了分層閾值函數,該函數結合了軟閾值和硬閾值函數去噪的優點,在保留了盡可能多的信號特征的情況下,使重構后的信號過渡更加平滑,其表達式為:
(9)
式中,σ為噪聲標準方差,j為小波分解層數,Tj為閾值。這里的N根據上文有關分解層數值的分析,取6為宜。這種方式減少了閾值處理的固定偏差,更好地保留了脈搏波信號在分解各層的特征,有效提升了降噪效果。
中值濾波是一種基于排序統計理論的能有效抑制噪聲,并對脈沖信號有良好濾除作用的非線性信號處理技術。這種濾波方式的基本原理是將數字圖像或數字序列中一點的值,用該點的一個鄰域中各點值的中值來代替[12],使周圍的像素值接近真實值,由此達到消除孤立噪聲點的效果。
中值濾波的數學表示為:

(10)

小波濾波在濾波處理上有獨特優勢,能對成分進行獨立分析,但不能對粗差的剔除做到完善。中值濾波在有效解決粗差的剔除上,具有獨特優勢[13]。因此,本研究采用了基于高頻中值濾波的小波濾波,用于處理采集的PPG信號,改善了僅采用小波濾波的局限性。
首先在改進算法中,對PPG信號用中值濾波處理,得到平滑且已經剔除高頻噪聲的信號。然后將處理后的信號進行小波分解,得到相對應的分解系數,再通過恰當的閾值函數[14]進行處理。最后將得到的系數重構,得到所需信號。
根據上述理論,可將基于高頻中值濾波的小波濾波流程進行歸納,見圖2。

圖2高頻中值濾波的小波濾波步驟
Fig.2The steps of wavelet filtering of high frequency median filter
具體步驟:
(1)對含有噪聲的信號作初步分析,由于信號中的噪聲主要以高頻系數形式出現,所以,使用中值濾波處理信號中的高頻系數部分;
(2)將得到的高頻系數與信號中的低頻系數加以重構,得到一個完成初步濾波的信號;
(3)在步驟(2)得到的信號上再使用小波濾波進行處理;
(4)選取合適的小波基,確定小波分解系數;
(5)對于步驟(4)中分解所得的分解系數,選擇合適的軟閾值或硬閾值進行閾值化處理,對各尺度小波系數去噪[15];
(6)對小波進行逆變換,將經軟硬閾值處理過的小波系數加以重構,最終得到基于中值濾波的小波濾波的特征信號。
為了驗證所提出方法的有效性,本研究利用表征觀察值與真值離散均方根誤差。其中,RMSE越小,信號處理結果越接近目標信號。
(11)
式中,N是樣本xn的個數。
同時利用所得的均方根誤差,計算信噪比,信噪比越大,說明混在信號里的噪聲越小。信噪比單位為分貝(dB)。

(12)
式中,Ps和Pn分別代表信號和噪聲的有效功率。
信號的頻譜圖表明了信號在不同頻率分量成分的大小,相比時域圖像提供了更豐富具體的頻域圖像。可以利用離散傅立葉變換實現時快速算法[16]即快速傅里葉變換(fast fourier transform, FFT),對所得信號進行頻譜分析,在要求精度不高的情況下,可在頻譜圖中讀出信號的成分,觀察頻譜圖可以看出高頻信號是否在經過中值濾波后被有效剔除。
(13)
式中,F(ω)為f(t)的像函數,f(t)為F(ω)的像原函數。
為了更好地比較不同濾波方式所達到的濾波效果,在使用從人體獲取的真實脈搏波進行數據處理之前,利用功能強大的MATLAB,使用仿真技術構建調試環境,對真實信號進行仿真[17-18]。通過仿真技術模擬實驗所獲得的參數結果,為真實數據的處理提供必要參考。
4.2.1信號仿真的步驟 本研究按照PPG信號的大致趨勢以及真實信號噪聲成分復雜的特征進行構建,使構建完畢的仿真信號與真實信號相比能夠獲得理想的均方根誤差值以及相關系數。
為了還原采樣時實測光電容積脈搏波信號在一個周期內的變化,在構建時,選擇了將正弦信號與余弦信號相疊加[19],以還原真實信號走勢,見圖3(a)。
由于PPG信號是心臟跳動所引發的波動信號,因此,該信號的主要頻率首先應根據心臟搏動所產生的心率來決定。臨床上,我們往往用脈率(pulse rate,PR)來替代,根據成人在靜息狀態PR參考值為60次/min,構建脈搏容積波頻率為1 Hz的仿真信號,見圖3(b)。
由于實測信號在測量中受到多種干擾因素影響,尤其是呼吸產生的干擾[20],進而使用高斯噪聲來還原這種干擾。高斯噪聲是一種較為泛用的噪聲模型,根據成人的呼吸頻率以及呼吸噪聲[21]的不規則性分布,在信號不同的區間上,分別疊加了均值為0,標準差分別為30、20、5、15、19、2的正態分布高斯噪聲,見圖3(c)。
疊加后,圖像呈現了較為明顯的不規則性。此外在測量過程中還存在工頻干擾以及測量設備傳輸過程中所產生的干擾,根據其特征,疊加了伽馬噪聲,該伽馬噪聲由三個服從指數分布的噪聲疊加而成,對此仿真圖像的構建具有十分顯著的效果,見圖3(d)。
仿真完成后,對構建完畢的仿真信號進行評價,將仿真信號與通過實測獲得的真實信號進行對比,經計算,得到RMSE為2.38,相關系數為0.928,符合實驗要求。

圖3PPG仿真信號疊加
Fig.3Superposition of simulation PPG signals
4.2.2特征信號的提取 在MATLAB中導入了構建的仿真信號。為將其干擾噪聲有效剔除,選擇了常見且實用的中值濾波。將仿真信號放入中值濾波器。中值濾波后的圖像見圖4。

圖4中值濾波效果圖
Fig.4Median filter rendering
從圖中可以看出,中值濾波對高頻噪聲有著較為顯著的濾波效果。對比原信號和濾波后的頻譜圖,在高頻區域內的噪聲已經被完全剔除,符合實驗目的。但是根據中值濾波后呈現出的頻譜圖像,見圖5,能夠發現低頻信號并未被有效過濾,如呼吸、心動等頻率低的干擾信號仍對特征信號的提取有一定影響。因此在中值濾波后再次選用了小波濾波。

圖5 原信號與中值濾波后的信號頻譜圖
利用小波濾波進行實驗時,首先需要選用合適的小波基和分解系數(N=6),之后對中值濾波后的信號進行分解,可以得到6層分解信號和位于第一層的中值濾波后信號,分解信號的圖像按照頻率由低到高依次排列。在閾值作用方式上,采用了能使圖像更平滑的軟閾值函數這一作用方式,同時也避免了硬閾值“一刀切”帶來的影響。在閾值的選取方面,使用了分層閾值函數進行閾值選取操作,能使濾波后信號的相似程度更高。根據式(11)計算可得,每一層閾值分別為99.6031、74.6282、62.7461、55.4213、50.3208、46.5034。這樣能使所得信號在相似度較高的基礎上,確保光滑度。
通過觀察中值濾波和組合濾波后信號的頻譜分析圖,可以看出,之前未被有效剔除的低頻干擾信號,在結合小波濾波后被有效去除,效果見圖6。綜上,此組合方法能有效濾除仿真光電容積脈搏波的干擾噪聲,可嘗試將其運用于真實的光電容積脈搏波信號濾波。

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

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

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