基于胸部生物阻抗技術開展對心臟血流動力學參數準確計算的前提是心阻抗微分信號特征點的準確識別。為提高心阻抗微分信號特征點識別的準確率，本文首先設計了自適應集合經驗模態和小波閾值降噪技術相結合的信號預處理方法，然后根據自適應集合經驗模態分解的結果，結合差分法和自適應分段技術定位了心阻抗微分信號中的 A、B、C 和 X 點。本研究以臨床上采集到的 30 例病理性心阻抗微分信號為例，對本文所提算法的特征檢測準確度進行檢驗。研究結果顯示，本文所提算法對信號特征的準確識別率總體可達 99.72%，進一步保證了基于胸部生物阻抗技術的心臟血流動力學參數的計算準確度。
引用本文： 張亞丹, 季忠, 譚霞, 申哲, 許戀姣. 用于無創心功能檢測的心阻抗微分信號處理. 生物醫學工程學雜志, 2019, 36(1): 50-58. doi: 10.7507/1001-5515.201804014 復制
引言
心阻抗微分信號（impedance cardiogram，ICG）是在心臟收縮、舒張過程中，反映胸部生物阻抗變化情況的信號，臨床上主要用于心功能檢查以得到心泵功能、心臟前后負荷及心肌收縮力等信息，其頻率范圍為 0.8～20 Hz[1]。對 ICG 信號的形態、幅值、波形寬度等特征進行分析，計算相關的血流動力學參數如心搏量、心排量、射血分數等，可預測或判斷心臟的功能性變化，而計算此類參數的前提是對左心室射血時間的準確獲取[2-4]。標準的 ICG 信號和心電信號（electrocardiogram，ECG）的特征對應關系如圖 1 所示[5]。z（t）為原始心阻抗信號，ICG 為 z（t）的微分，其中 ICG 信號中，A 點對應于心室收縮的起點；B 點為主動脈瓣開放的點，即心室射血起點；X 點是主動脈瓣關閉的點，即心室射血終點；左心室射血時間即為 B 點到 X 點的時間間隔；O 點為二尖瓣打開點。而 ECG 信號中，QRS 波反映心臟受到電刺激后，心室去極化的過程；而 T 波反映心室復極化的過程。ECG 信號的 Q 點到 ICG 信號的 B 點的時間為心臟的射血前期；ECG 信號的 Q 波起點到 ICG 信號的峰點 C 的時間為左室功能指數[6-7]。
圖1
ICG 信號和 ECG 信號特征關系圖
Figure1.
The corresponding features between ECG and ICG signals
基于 ICG 信號和 ECG 信號分析心臟血流動力學參數的方法具有連續、無創、價格低廉、操作簡便等優勢，但現階段該方法并未在臨床上得以廣泛應用，這是由于 ICG 信號的微弱性，且在病理或受干擾狀態下的復雜性導致其特征點較難定位，繼而嚴重影響相關心臟血流動力學參數的準確性及可靠性的緣故[8-9]，因此如何保證 ICG 信號特征的準確識別，特別是病理狀態下的 ICG 信號的特征識別是現階段亟待解決的問題。
目前，對 ICG 信號的預處理方法主要有基于時域或頻域的濾波法，以及基于時頻域的小波法、集合經驗模態分解法（ensemble empirical mode decomposition，EEMD）等。Hu 等[1]提出一種基于最小均方誤差的自適應濾波算法來濾除 ICG 信號中由于呼吸運動引起的漂移噪聲。Janu?auskas 等[10]提出了利用 EEMD 的方法對 ICG 信號進行降噪，通過設計實驗證明了 EEMD 法對 ICG 信號的降噪效果要優于離散小波變換法和線性濾波法。Chabchoub 等[11]利用離散小波變換的方法對 ICG 信號進行降噪，通過和基于最小二乘的卷積擬合濾波法（Savitzky-Golay 濾波法）、基于均值濾波的降噪算法進行比較，證明了小波降噪算法性能最好。Choudhari 等[12]針對 ICG 信號中典型的基線漂移噪聲，分別利用小波包分解和卡爾曼濾波的方法對 ICG 信號進行降噪，結果證明，用 db4 小波對信號進行小波包分解后的降噪效果更佳。劉珊[13]選用 bior3.5 小波作為小波基函數對 ICG 信號進行小波分解，針對低頻噪聲所在的階層，對其近似系數求平均值以達到降噪的目的；對于高頻噪聲分布的階層，對其低頻系數和高頻系數分別進行閾值量化處理，最終重構各級處理后的小波系數獲取降噪后信號。單一的時域分析或頻域分析的方法不能較好地濾除頻率成分較復雜的噪聲，而時頻域分析法，如小波法和 EEMD 法等可對信號進行多分辨率分解，從而可實現對頻帶范圍重疊于有用信號頻帶范圍的較復雜的噪聲成分的濾除。但是以上的時頻域分析法都不可避免地存在主觀經驗誤差，如：小波法中小波基的選擇、EEMD 法中加入的白噪聲量和集合平均次數的確定等，都需要參考主觀經驗。信號降噪效果的好壞會直接影響到后續信號各特征點的準確識別，因此如何設計出一套有效的信號降噪算法至關重要。
對于 ICG 信號的特征點，即 A 點、B 點、C 點和 X 點的識別方法可以分為直接法和間接法兩類。直接法即通過分析 ICG 信號自身的特點，對其特征信息進行識別；間接法為結合同步采集的 ECG 信號，參考 ECG 信號和 ICG 信號特征點的對應關系，利用 ECG 信號的特征點標校 ICG 信號的特征點[5]。Podtaev 等[7]提出了一種小波差分算法完成了對 ICG 信號的特征點的檢測；趙越[14]首先利用求最大值的方法求出 C 點，利用 B 點落后于 ECG 信號 R 點和 B 點與 C 點的幅值關系定位 B 點，以 T 波終點標校 X 點，最終完成 ICG 信號特征點的定位；趙云冬等[6]用 bior3.7 小波對 ICG 信號進行降噪和特征點檢測，對健康人和心血管疾病患者（以下簡稱為患者）的 ICG 信號進行預處理，實現了 C 點、B 點和 X 點等特征點的識別，指出差分閾值法對患者信號的檢測準確率為 94.38%～95.49%，對健康人的準確檢出率為 96.37%～97.65%，小波變換法對患者的檢測準確率為 97.13%～98.10%，對健康人的準確檢出率為 97.58%～98.51%。
綜上所述，已有的很多研究對 ICG 信號的降噪和特征點提取都傾向于采用小波分解的方式，但是小波分解方式中，在進行小波基的選擇時，除需考量小波基自身的特征外，還需考慮小波基與待處理信號的相似性等因素，因此會不可避免地引入主觀經驗誤差，且根據現有相關研究，小波變換法對患者 ICG 信號的檢測準確率最高僅為 98.10%，該值仍有很大的上升空間。因此，本文提出了以自適應集合經驗模態分解（adaptive ensemble empirical mode decomposition，AEEMD）算法為基礎的 ICG 信號降噪和特征點檢測算法。AEEMD 算法是由經驗模態分解法（empirical mode decomposition，EMD）發展起來的，同小波變換一樣，AEEMD 算法也具有多分辨率的特點，且避免了小波基選擇的難題[15]。本文針對 ICG 信號和噪聲信號各自的頻率分布特點，在對信號進行 AEEMD 分解的基礎上，分別設計不同的方法去除不同頻率成分的噪聲；根據 AEEMD 算法的分解結果，選擇特定階層的信號形成特征檢測層，依據特征檢測層和降噪后信號間各特征點的對應位置關系，定位出最終的 A、B、C 和 X 點，以期提高特征點識別的準確率，從而確保相關血流動力學參數的準確計算，達到改善基于胸部生物阻抗技術的無創心功能檢測分析儀的魯棒性的目的，為醫生提供更加可靠的診斷依據。
1 原理和方法
1.1 EMD 和 EEMD 原理
EMD 是一種無需先驗知識，無需基函數的可以以自身信號特點為依托，把復合頻率的非平穩信號自適應地分解成從高頻至低頻的一系列固有模態函數（intrinsic mode function，IMF）的時頻分析方法[16]。
對信號進行 EMD 分解時，首先應找出信號的極值點，分別由極大值點和極小值點確定出該信號的上下包絡線，然后對其相加求均值，用原信號減去該均值即得第一階 IMF（以符號 IMF1 表示）；為避免求包絡過程中由于信號局部某處凹凸引起錯檢極值點的誤差，并使求得的波形更加對稱，需對檢出的 IMF1 進行檢驗，判斷其是否為理想狀況下的 IMF；若是，則繼續確定第二階 IMF（以符號 IMF2 表示），以此類推，否則將第一次確定的理想狀態下的 IMF1 作為原始信號，重復上述求極值、求包絡、求均值和相減的操作，直至得到滿足 IMF 條件的 IMF1 為止。用原始信號減去 IMF1，并將其作為原始信號重復以上操作，以求得 IMF2、IMF3，
，求每階 IMF 時，對應原始信號與該階 IMF 求差可得到對應階層的剩余信號，若剩余信號（以符號 r 表示）足夠小或剩余信號無極值存在時，終止分解過程。
EEMD 是在 EMD 的基礎上發展起來的，在該方法中，為保證待分析信號的時域連續性，首先向待分析信號中疊加不同數量的白噪聲，以達到避免插值過程中過沖和欠沖發生的目的；然后，對各個疊加了白噪聲之后的信號進行 EMD 分解，分別得到對應的各階 IMF1、IMF2、IMF3，
；最后，對各對應階層的 IMF 進行疊加平均，從而得到最終的 IMF1、IMF2、IMF3，
。該方法有效地抑制了 EMD 中的模態混疊現象。在 EEMD 方法中，加入的白噪聲量和集合平均次數一般由經驗確定，這就會不可避免地引入一些主觀經驗誤差[17]。
1.2 AEEMD 原理
AEEMD 算法依據信號本身的特點自適應地確定 EEMD 中所需的兩個參數：待加入的白噪聲量和集合平均次數。這樣處理不僅解決了 EMD 中存在的模態混疊問題，也避免了 EEMD 中存在的依據經驗設置兩參數而存在的經驗誤差問題。該方法對信號進行分解的具體步驟如下：
（1）對待處理信號（以符號 x（n）表示）進行 EMD 分解，得到一系列頻率范圍與分解階數負相關的各階 IMF[18]。
（2）取前兩階 IMF 進行疊加作為信號高頻噪聲的估計，其幅值標準差（standard deviation，SD）以符號 SD_HF 表示；原始信號的幅值標準差以符號 SD_X 表示。
（3）確定待加入的白噪聲的幅值標準差（以符號 SD_WD 表示）和集合平均次數（以符號 N 表示）。合適的白噪聲量應使信號低頻部分的極值點分布均勻，對高頻部分極值點的分布幾乎不產生影響，從而保證分解結果的可靠性。基于概率統計的知識，添加的 SD_WD 和 N 分別如式（1）、（2）所示[19]：
 |
 |
其中，SD_WD 取 SD_HF 的 1/3，ε 為信號分解的可容許誤差，一般 ε 為 0.01 即可滿足要求。
（4）對信號進行 EEMD 分解。
1.3 ICG 信號的預處理和特征點識別算法
針對 ICG 信號在采集過程中受到的不同類型的噪聲干擾，為盡量避免信號多尺度分解時造成的誤差，本文選用 AEEMD 算法對信號進行多分辨率分解。在保證噪聲被有效濾除的情況下，若對分解后的高頻信號直接進行取舍操作或簡單地閾值處理，有用信號的能量不能較好地保存下來[20]，故本文對分解后的各級信號中的高頻部分采用小波閾值法（wavelet thresholding，WT）進行二級降噪。WT 法不僅可對信號進行多分辨率分解，而且在計算效率上較優于 AEEMD 算法。此外，在本文中，僅對 AEEMD 算法分解后的高頻部分做小波分解時，固定小波基帶來的誤差在整體信號的降噪中可忽略不計。因此，本文提出了 AEEMD 算法和 WT 算法相結合法（AEEMD and WT，AEEMD_WT）對信號進行預處理。
1.3.1 基于 AEEMD_WT 的信號降噪
采用 AEEMD_WT 法對信號降噪的過程中，首先，對信號進行 AEEMD 分解；然后，根據分解結果中各階 IMF 的頻率范圍，選擇固定階數的 IMF 進行疊加，如本文信號的采樣頻率為 1 000 Hz，取前四階 IMF 進行疊加，根據各階 IMF 頻率按近似 2 的負冪次方遞減的規律[13]，疊加后的前四階信號的頻率范圍主要集中在 62.5～1 000 Hz，為更大程度地保留前四階 IMF 中的有用信號，本文結合 WT 法對該疊加層信號中較高頻噪聲進行濾除。WT 算法降噪法中，選用 coif4 作為小波基函數，分解層數為 6 層，閾值采用伯吉馬薩特（Brige-Massart）策略確定，用軟閾值法對各級小波系數進行處理，最終重構各級系數得到降噪后的前四階 IMF 疊加信號（以符號 IMFS4 表示）。IMF9 及其以后各階 IMF 主要為信號中所含較低頻噪聲所在的頻段，對其直接舍去；最后，將 IMFS4 和 IMF5、IMF6、IMF7、IMF8 進行疊加，重構降噪后的信號（以符號 c（n）表示）。
1.3.2 基于 AEEMD、差分法和自適應分段技術的 ICG 信號特征點識別
考慮到 ICG 信號的微弱性及病理狀態下 ICG 信號的多樣性，本文提出了一種結合 AEEMD、差分法和自適應分段技術（AEEMD，difference and adaptive segmentation，AEEMD_D_AS）的方法對 ICG 信號進行特征點識別。考慮到本文信號 1 000 Hz 的采樣頻率，為保證對信號 AEEMD 分解的效率及分解結果的可靠性，本文以 7 000 個點為一個處理單位，對 ICG 信號進行特征點的檢測，步驟如下：
（1）原始信號 x（n）進行 AEEMD 分解的結果，如式（3）所示：
 |
其中，IMFm，1，IMFm+1，2，
，IMFm+n-1，n 為信號頻率占主要成分的 IMF 階層。
（2）為避免干擾波形的影響，突出待檢特征點，本文選取各特征波占主要頻率成分的 IMF 層進行合成，作為 A、B、C、X 點的檢測層信號（以符號 d（n）表示）。
 |
（3）為避免各節拍信號特征的漏檢、誤檢，本文對信號進行分段處理，以保證各段信號中僅包含一個節拍長度的信號（分段長度以符號 len 表示），且設計反饋機制：考慮到病理狀態下，各節拍的長度可能發生變化，本文以檢測到的連續三個波形周期的方差之和取極小值或該方差取到一定范圍內的最小值為終止條件，自適應地更新 len，len 的每次變化量（以符號 sp 表示）為 0.1 倍的采樣頻率。參考 ECG 信號的周期及 ECG 信號和 ICG 信號的對應關系，len 的變化范圍設計為 0.4～1.3 倍的采樣頻率。
（4）對 d（n）求差分（以符號 d_d（n）表示該差分信號），以增強待檢波形的能量占比。
（5）初始化 n = 0。
（6）在 d_d ［（1 + n·sp）：（n + 1）·len］信號中，計算最小值（符號記為：min_value），且記錄該最小值所在位置（符號記為：min_index），并據此定位最小值左側第一個最大值（符號記為：max_value），且記錄其所在位置（符號記為：max_index）。
（7）經反復實驗，將步驟（6）中定位出的極值對中間的過零點作為準 C 點對應位置（符號記為：QC_index）。為避免噪聲的干擾，本文設計了 C 點的校準程序。在降噪后的信號 c［QC_index –（QC_index – max_index）/2：QC_index + （min_index – QC_index）/2］上檢測最值點，該最值點若為準 C 點所在位置，則準 C 點位置即為待識別的 C 點所在位置；若非準 C 點所在位置，考慮到檢測層信號并不能完全反映信號特征點的位置信息，故標記降噪后信號中檢測到的最值點為最終的 C 點，分別記錄 C 點所在位置（符號記為：C_index）及幅值（符號記為：xC）。
（8）為提高檢測效率，在 C 點左側一定范圍內檢測 A 點，考慮到波形的對稱性，確定 A 點的檢測范圍為 d_d{［（2·max_index – QC_index –（QC_index – max_index）/2］：［2·max_index – QC_index + （QC_index – max_index）/2）］}，在該范圍內檢測距離 C 點最近的過零點位置（符號記為：QA_index），該位置即為準 A 點所在位置。
（9）判斷點 c［QA_index，x（QA_index）］是否為最值點，若是，則上步中所求準 A 點即為最終 A 點；否則，在 c［QA_index –（max_index – QA_index）/2：QA_index +（max_index – QA_index）/2］范圍內檢測最值點，該最值點所在位置即為最終的 A 點所在位置（符號記為 A_index）。
（10）最終的 B 點在 c（A_index：max_index）中幅值等于 0.15·xC 位置處。
（11）檢測 X 點。在 d_d［（1 + n·sp）：（n + 1）· sp］上，C 點右側一定范圍內檢測 X 點。考慮到算法效率，縮小 X 點的檢測范圍為 d_d{［2·min_index – QC_index –（min_index – QC_index）/2］：［2·min_index – QC_index +（min_index – QC_index）/2）］}，在該范圍內檢測距離 min_index 最近的一個過零點，對應位置為準 X 點所在位置（符號記為：QX_index）。
（12）判斷準 X 點 c［QX_index，x（QX_index）］是否為極值點，若是，則此點即為 X 點；否則，在 c［QX_index –（QX_index – min_index）/2：QX_index + （QX_index – min_index）/2］范圍內檢測極值點，將該極值點作為最終的 X 點。
（13）n = n + 1。重復步驟 （1）～（11），直到檢完整段信號。
2 實驗及結果
2.1 實驗數據來源
本次實驗中所用數據由重慶市武警重慶總隊醫院的心功能科室提供，數據的采集設備為心功能無創檢測分析儀（KF_ICG-303 型，重慶科發醫療器械有限公司，重慶）[21]。本文選取 30 例形態差異較明顯的患者的信號，對本文所提算法的有效性進行評估。其中，30 例患者的年齡分布于 31～88 歲，男性患者 17 例，女性患者 13 例。針對各例患者，分別選取時長約 5 min 的信號，對其特征信息進行識別、統計。
2.2 實驗結果及分析
本文隨機選取一例患者信號（采樣頻率為 1 000 Hz），對實驗過程及結果進行詳細論述如下。
部分原始信號如圖 2 所示。
圖2
部分原始信號
Figure2.
Unprocessed signal fraction
該信號夾雜噪聲，會影響對信號的特征提取結果，故首先對信號進行降噪處理。對信號進行 AEEMD 分解，AEEMD 分解后的各階 IMF 如圖 3 所示。
圖3
待處理信號 AEEMD 分解結果
Figure3.
Decomposed results in AEEMD of unprocessed signal
選取較高頻噪聲占主要成分的 IMFs4 信號進行疊加，并用 WT 法對其進行降噪處理，對基線漂移占主要成分的 IMF9、IMF10、IMF11 直接舍去，重構 WT 算法降噪后的 IMFs4 信號及有用信號占主要成分的 IMF5、IMF6、IMF7、IMF8 信號，最終對 ICG 信號的降噪結果如圖 4 所示。
圖4
降噪后信號
Figure4.
The de-noised signal
濾波前后，信號的頻譜變化如圖 5 所示。
圖5
信號降噪前后信號頻譜變化
Figure5.
The frequency spectrum of signal before and after de-noise
整體來看，降噪后信號中較高頻噪聲均被有效濾除，根據降噪前后信號的 20 Hz 以內的頻譜分布情況，較低頻噪聲和部分重疊頻率噪聲對應位置的頻率也得以明顯衰減，證明了 AEEMD_WT 法對 ICG 信號中多種類型的噪聲濾除的有效性。
此外，經反復實驗論證，在 IMF6、IMF7 疊加完進行差分后的信號中，各節拍信號的過零點位置即為對應 ICG 信號中的 C 點位置，故將此作為特征波檢測層。特征檢測層和降噪后信號的特征對應關系如圖 6 所示。
圖6
檢測層和降噪后信號特征對應關系
Figure6.
The corresponding relationship between the de-noised signal and the detection signal
如圖 6 所示，檢測層差分信號 d_d（n）中的橫虛線代表該信號的過零點位置，豎虛線表示了 d_d（n）和降噪后的 ICG 信號 c（n）的特征點的對應關系。在單個分段長度內，c（n）中 C 點位于 d_d（n）中的最大值（極大值）和最小值（極小值）間過零點的位置，本文選擇在 c（n）上的該極值對范圍內檢測極值點，以此確定 C_index；在 d_d（n）信號中，C_index 左側第一個過零點位置應為理想情況下的 A_index，C_index 右側第一個過零點位置為理想情況下 X_index，為避免噪聲的干擾，在 c（n）中，以 d_d（n）信號中的極大值、極小值位置為參考，確立 A 點和 X 點的檢測范圍，檢測對應的極值點以避免誤檢。特征點檢測效果如圖 7 所示。
圖7
特征識別結果
Figure7.
The result of characteristic recognition
參照標準 ICG 信號各特征點的位置關系，并經過和有多年相關臨床經驗的武警重慶總隊醫院的心功能科室醫生反復討論分析，上述 ICG 信號的各特征點均可認為被有效檢出。研究結果顯示，采用本文所提算法對臨床上大量不同病癥類型的信號進行分析處理，均取得了滿意的結果。最終的特征點檢測統計結果如表 1 所示。需要指出的是，本文所提出的特征點識別算法中，A 點、B 點、X 點的檢測均有賴于 C 點的成功檢出，即任一節拍中 C 點出現漏檢或誤檢，則整個節拍都將漏檢或誤檢；所謂誤檢，為錯檢和多檢的節拍總數。考慮到各特征信息的位置關系及基于此的算法設計原理，統計時，本文選擇以整個節拍為單位（任一特征點出現誤檢，即視為此節拍誤檢），參照公式（5）、（6）、（7），分別計算漏檢率、誤檢率和準確率，最終的統計結果如表 1 所示。
表1
特征點檢測結果統計
Table1.
The statistical result of feature points detected
 |
 |
 |
特征檢測過程中，多處漏檢和誤檢均發生在信號自身特征缺失處。分析原因可能是由于 ICG 信號的微弱性使其在采集過程中極易受到噪聲的干擾，特別是噪聲頻率重疊于有用信號頻率的情況，導致信號特征被淹沒；另一方面，受試者在測量過程中偶見的保護性反射行為，如咳嗽等，均會影響信號形態特征；此外，受試者均有不同類型的心血管疾病，疾病的嚴重程度和類型會影響 ICG 信號的形態特征，如幅值驟減、信號變形、甚至部分特征點消失等。
3 結論
臨床上，基于胸部生物阻抗技術的心臟血流動力學分析的可靠性嚴重受限于 ICG 信號特征的準確識別，為保證心功能無創檢測儀更好地發揮其臨床應用價值，本文首先考慮到 ICG 信號所含噪聲的多樣性，提出了一種 AEEMD_WT 降噪算法，成功濾除了信號中較高頻噪聲、低頻噪聲及肌電干擾、電極移動等造成的和有用信號頻率有重疊的噪聲，為特征點的準確識別奠定了基礎；繼而，考慮到病理性 ICG 信號的形態多變性，設計了可自適應更新檢測節拍長度的 AEEMD_D_AS 特征檢測算法完成了對 A 點、B 點、C 點和 X 點的檢測；最后，以臨床上隨機抽樣到的 30 例病理性 ICG 信號對降噪及特征識別算法進行驗證，最終檢測的準確率高達 99.72%，進一步保證了基于胸部生物阻抗技術的心臟血流動力學參數的計算準確度。
引言
心阻抗微分信號（impedance cardiogram，ICG）是在心臟收縮、舒張過程中，反映胸部生物阻抗變化情況的信號，臨床上主要用于心功能檢查以得到心泵功能、心臟前后負荷及心肌收縮力等信息，其頻率范圍為 0.8～20 Hz[1]。對 ICG 信號的形態、幅值、波形寬度等特征進行分析，計算相關的血流動力學參數如心搏量、心排量、射血分數等，可預測或判斷心臟的功能性變化，而計算此類參數的前提是對左心室射血時間的準確獲取[2-4]。標準的 ICG 信號和心電信號（electrocardiogram，ECG）的特征對應關系如圖 1 所示[5]。z（t）為原始心阻抗信號，ICG 為 z（t）的微分，其中 ICG 信號中，A 點對應于心室收縮的起點；B 點為主動脈瓣開放的點，即心室射血起點；X 點是主動脈瓣關閉的點，即心室射血終點；左心室射血時間即為 B 點到 X 點的時間間隔；O 點為二尖瓣打開點。而 ECG 信號中，QRS 波反映心臟受到電刺激后，心室去極化的過程；而 T 波反映心室復極化的過程。ECG 信號的 Q 點到 ICG 信號的 B 點的時間為心臟的射血前期；ECG 信號的 Q 波起點到 ICG 信號的峰點 C 的時間為左室功能指數[6-7]。
圖1
ICG 信號和 ECG 信號特征關系圖
Figure1.
The corresponding features between ECG and ICG signals
基于 ICG 信號和 ECG 信號分析心臟血流動力學參數的方法具有連續、無創、價格低廉、操作簡便等優勢，但現階段該方法并未在臨床上得以廣泛應用，這是由于 ICG 信號的微弱性，且在病理或受干擾狀態下的復雜性導致其特征點較難定位，繼而嚴重影響相關心臟血流動力學參數的準確性及可靠性的緣故[8-9]，因此如何保證 ICG 信號特征的準確識別，特別是病理狀態下的 ICG 信號的特征識別是現階段亟待解決的問題。
目前，對 ICG 信號的預處理方法主要有基于時域或頻域的濾波法，以及基于時頻域的小波法、集合經驗模態分解法（ensemble empirical mode decomposition，EEMD）等。Hu 等[1]提出一種基于最小均方誤差的自適應濾波算法來濾除 ICG 信號中由于呼吸運動引起的漂移噪聲。Janu?auskas 等[10]提出了利用 EEMD 的方法對 ICG 信號進行降噪，通過設計實驗證明了 EEMD 法對 ICG 信號的降噪效果要優于離散小波變換法和線性濾波法。Chabchoub 等[11]利用離散小波變換的方法對 ICG 信號進行降噪，通過和基于最小二乘的卷積擬合濾波法（Savitzky-Golay 濾波法）、基于均值濾波的降噪算法進行比較，證明了小波降噪算法性能最好。Choudhari 等[12]針對 ICG 信號中典型的基線漂移噪聲，分別利用小波包分解和卡爾曼濾波的方法對 ICG 信號進行降噪，結果證明，用 db4 小波對信號進行小波包分解后的降噪效果更佳。劉珊[13]選用 bior3.5 小波作為小波基函數對 ICG 信號進行小波分解，針對低頻噪聲所在的階層，對其近似系數求平均值以達到降噪的目的；對于高頻噪聲分布的階層，對其低頻系數和高頻系數分別進行閾值量化處理，最終重構各級處理后的小波系數獲取降噪后信號。單一的時域分析或頻域分析的方法不能較好地濾除頻率成分較復雜的噪聲，而時頻域分析法，如小波法和 EEMD 法等可對信號進行多分辨率分解，從而可實現對頻帶范圍重疊于有用信號頻帶范圍的較復雜的噪聲成分的濾除。但是以上的時頻域分析法都不可避免地存在主觀經驗誤差，如：小波法中小波基的選擇、EEMD 法中加入的白噪聲量和集合平均次數的確定等，都需要參考主觀經驗。信號降噪效果的好壞會直接影響到后續信號各特征點的準確識別，因此如何設計出一套有效的信號降噪算法至關重要。
對于 ICG 信號的特征點，即 A 點、B 點、C 點和 X 點的識別方法可以分為直接法和間接法兩類。直接法即通過分析 ICG 信號自身的特點，對其特征信息進行識別；間接法為結合同步采集的 ECG 信號，參考 ECG 信號和 ICG 信號特征點的對應關系，利用 ECG 信號的特征點標校 ICG 信號的特征點[5]。Podtaev 等[7]提出了一種小波差分算法完成了對 ICG 信號的特征點的檢測；趙越[14]首先利用求最大值的方法求出 C 點，利用 B 點落后于 ECG 信號 R 點和 B 點與 C 點的幅值關系定位 B 點，以 T 波終點標校 X 點，最終完成 ICG 信號特征點的定位；趙云冬等[6]用 bior3.7 小波對 ICG 信號進行降噪和特征點檢測，對健康人和心血管疾病患者（以下簡稱為患者）的 ICG 信號進行預處理，實現了 C 點、B 點和 X 點等特征點的識別，指出差分閾值法對患者信號的檢測準確率為 94.38%～95.49%，對健康人的準確檢出率為 96.37%～97.65%，小波變換法對患者的檢測準確率為 97.13%～98.10%，對健康人的準確檢出率為 97.58%～98.51%。
綜上所述，已有的很多研究對 ICG 信號的降噪和特征點提取都傾向于采用小波分解的方式，但是小波分解方式中，在進行小波基的選擇時，除需考量小波基自身的特征外，還需考慮小波基與待處理信號的相似性等因素，因此會不可避免地引入主觀經驗誤差，且根據現有相關研究，小波變換法對患者 ICG 信號的檢測準確率最高僅為 98.10%，該值仍有很大的上升空間。因此，本文提出了以自適應集合經驗模態分解（adaptive ensemble empirical mode decomposition，AEEMD）算法為基礎的 ICG 信號降噪和特征點檢測算法。AEEMD 算法是由經驗模態分解法（empirical mode decomposition，EMD）發展起來的，同小波變換一樣，AEEMD 算法也具有多分辨率的特點，且避免了小波基選擇的難題[15]。本文針對 ICG 信號和噪聲信號各自的頻率分布特點，在對信號進行 AEEMD 分解的基礎上，分別設計不同的方法去除不同頻率成分的噪聲；根據 AEEMD 算法的分解結果，選擇特定階層的信號形成特征檢測層，依據特征檢測層和降噪后信號間各特征點的對應位置關系，定位出最終的 A、B、C 和 X 點，以期提高特征點識別的準確率，從而確保相關血流動力學參數的準確計算，達到改善基于胸部生物阻抗技術的無創心功能檢測分析儀的魯棒性的目的，為醫生提供更加可靠的診斷依據。
1 原理和方法
1.1 EMD 和 EEMD 原理
EMD 是一種無需先驗知識，無需基函數的可以以自身信號特點為依托，把復合頻率的非平穩信號自適應地分解成從高頻至低頻的一系列固有模態函數（intrinsic mode function，IMF）的時頻分析方法[16]。
對信號進行 EMD 分解時，首先應找出信號的極值點，分別由極大值點和極小值點確定出該信號的上下包絡線，然后對其相加求均值，用原信號減去該均值即得第一階 IMF（以符號 IMF1 表示）；為避免求包絡過程中由于信號局部某處凹凸引起錯檢極值點的誤差，并使求得的波形更加對稱，需對檢出的 IMF1 進行檢驗，判斷其是否為理想狀況下的 IMF；若是，則繼續確定第二階 IMF（以符號 IMF2 表示），以此類推，否則將第一次確定的理想狀態下的 IMF1 作為原始信號，重復上述求極值、求包絡、求均值和相減的操作，直至得到滿足 IMF 條件的 IMF1 為止。用原始信號減去 IMF1，并將其作為原始信號重復以上操作，以求得 IMF2、IMF3，，求每階 IMF 時，對應原始信號與該階 IMF 求差可得到對應階層的剩余信號，若剩余信號（以符號 r 表示）足夠小或剩余信號無極值存在時，終止分解過程。
EEMD 是在 EMD 的基礎上發展起來的，在該方法中，為保證待分析信號的時域連續性，首先向待分析信號中疊加不同數量的白噪聲，以達到避免插值過程中過沖和欠沖發生的目的；然后，對各個疊加了白噪聲之后的信號進行 EMD 分解，分別得到對應的各階 IMF1、IMF2、IMF3，；最后，對各對應階層的 IMF 進行疊加平均，從而得到最終的 IMF1、IMF2、IMF3，。該方法有效地抑制了 EMD 中的模態混疊現象。在 EEMD 方法中，加入的白噪聲量和集合平均次數一般由經驗確定，這就會不可避免地引入一些主觀經驗誤差[17]。
1.2 AEEMD 原理
AEEMD 算法依據信號本身的特點自適應地確定 EEMD 中所需的兩個參數：待加入的白噪聲量和集合平均次數。這樣處理不僅解決了 EMD 中存在的模態混疊問題，也避免了 EEMD 中存在的依據經驗設置兩參數而存在的經驗誤差問題。該方法對信號進行分解的具體步驟如下：
（1）對待處理信號（以符號 x（n）表示）進行 EMD 分解，得到一系列頻率范圍與分解階數負相關的各階 IMF[18]。
（2）取前兩階 IMF 進行疊加作為信號高頻噪聲的估計，其幅值標準差（standard deviation，SD）以符號 SD_HF 表示；原始信號的幅值標準差以符號 SD_X 表示。
（3）確定待加入的白噪聲的幅值標準差（以符號 SD_WD 表示）和集合平均次數（以符號 N 表示）。合適的白噪聲量應使信號低頻部分的極值點分布均勻，對高頻部分極值點的分布幾乎不產生影響，從而保證分解結果的可靠性。基于概率統計的知識，添加的 SD_WD 和 N 分別如式（1）、（2）所示[19]：
|
|
其中，SD_WD 取 SD_HF 的 1/3，ε 為信號分解的可容許誤差，一般 ε 為 0.01 即可滿足要求。
（4）對信號進行 EEMD 分解。
1.3 ICG 信號的預處理和特征點識別算法
針對 ICG 信號在采集過程中受到的不同類型的噪聲干擾，為盡量避免信號多尺度分解時造成的誤差，本文選用 AEEMD 算法對信號進行多分辨率分解。在保證噪聲被有效濾除的情況下，若對分解后的高頻信號直接進行取舍操作或簡單地閾值處理，有用信號的能量不能較好地保存下來[20]，故本文對分解后的各級信號中的高頻部分采用小波閾值法（wavelet thresholding，WT）進行二級降噪。WT 法不僅可對信號進行多分辨率分解，而且在計算效率上較優于 AEEMD 算法。此外，在本文中，僅對 AEEMD 算法分解后的高頻部分做小波分解時，固定小波基帶來的誤差在整體信號的降噪中可忽略不計。因此，本文提出了 AEEMD 算法和 WT 算法相結合法（AEEMD and WT，AEEMD_WT）對信號進行預處理。
1.3.1 基于 AEEMD_WT 的信號降噪
采用 AEEMD_WT 法對信號降噪的過程中，首先，對信號進行 AEEMD 分解；然后，根據分解結果中各階 IMF 的頻率范圍，選擇固定階數的 IMF 進行疊加，如本文信號的采樣頻率為 1 000 Hz，取前四階 IMF 進行疊加，根據各階 IMF 頻率按近似 2 的負冪次方遞減的規律[13]，疊加后的前四階信號的頻率范圍主要集中在 62.5～1 000 Hz，為更大程度地保留前四階 IMF 中的有用信號，本文結合 WT 法對該疊加層信號中較高頻噪聲進行濾除。WT 算法降噪法中，選用 coif4 作為小波基函數，分解層數為 6 層，閾值采用伯吉馬薩特（Brige-Massart）策略確定，用軟閾值法對各級小波系數進行處理，最終重構各級系數得到降噪后的前四階 IMF 疊加信號（以符號 IMFS4 表示）。IMF9 及其以后各階 IMF 主要為信號中所含較低頻噪聲所在的頻段，對其直接舍去；最后，將 IMFS4 和 IMF5、IMF6、IMF7、IMF8 進行疊加，重構降噪后的信號（以符號 c（n）表示）。
1.3.2 基于 AEEMD、差分法和自適應分段技術的 ICG 信號特征點識別
考慮到 ICG 信號的微弱性及病理狀態下 ICG 信號的多樣性，本文提出了一種結合 AEEMD、差分法和自適應分段技術（AEEMD，difference and adaptive segmentation，AEEMD_D_AS）的方法對 ICG 信號進行特征點識別。考慮到本文信號 1 000 Hz 的采樣頻率，為保證對信號 AEEMD 分解的效率及分解結果的可靠性，本文以 7 000 個點為一個處理單位，對 ICG 信號進行特征點的檢測，步驟如下：
（1）原始信號 x（n）進行 AEEMD 分解的結果，如式（3）所示：
|
其中，IMFm，1，IMFm+1，2，，IMFm+n-1，n 為信號頻率占主要成分的 IMF 階層。
（2）為避免干擾波形的影響，突出待檢特征點，本文選取各特征波占主要頻率成分的 IMF 層進行合成，作為 A、B、C、X 點的檢測層信號（以符號 d（n）表示）。
|
（3）為避免各節拍信號特征的漏檢、誤檢，本文對信號進行分段處理，以保證各段信號中僅包含一個節拍長度的信號（分段長度以符號 len 表示），且設計反饋機制：考慮到病理狀態下，各節拍的長度可能發生變化，本文以檢測到的連續三個波形周期的方差之和取極小值或該方差取到一定范圍內的最小值為終止條件，自適應地更新 len，len 的每次變化量（以符號 sp 表示）為 0.1 倍的采樣頻率。參考 ECG 信號的周期及 ECG 信號和 ICG 信號的對應關系，len 的變化范圍設計為 0.4～1.3 倍的采樣頻率。
（4）對 d（n）求差分（以符號 d_d（n）表示該差分信號），以增強待檢波形的能量占比。
（5）初始化 n = 0。
（6）在 d_d ［（1 + n·sp）：（n + 1）·len］信號中，計算最小值（符號記為：min_value），且記錄該最小值所在位置（符號記為：min_index），并據此定位最小值左側第一個最大值（符號記為：max_value），且記錄其所在位置（符號記為：max_index）。
（7）經反復實驗，將步驟（6）中定位出的極值對中間的過零點作為準 C 點對應位置（符號記為：QC_index）。為避免噪聲的干擾，本文設計了 C 點的校準程序。在降噪后的信號 c［QC_index –（QC_index – max_index）/2：QC_index + （min_index – QC_index）/2］上檢測最值點，該最值點若為準 C 點所在位置，則準 C 點位置即為待識別的 C 點所在位置；若非準 C 點所在位置，考慮到檢測層信號并不能完全反映信號特征點的位置信息，故標記降噪后信號中檢測到的最值點為最終的 C 點，分別記錄 C 點所在位置（符號記為：C_index）及幅值（符號記為：xC）。
（8）為提高檢測效率，在 C 點左側一定范圍內檢測 A 點，考慮到波形的對稱性，確定 A 點的檢測范圍為 d_d{［（2·max_index – QC_index –（QC_index – max_index）/2］：［2·max_index – QC_index + （QC_index – max_index）/2）］}，在該范圍內檢測距離 C 點最近的過零點位置（符號記為：QA_index），該位置即為準 A 點所在位置。
（9）判斷點 c［QA_index，x（QA_index）］是否為最值點，若是，則上步中所求準 A 點即為最終 A 點；否則，在 c［QA_index –（max_index – QA_index）/2：QA_index +（max_index – QA_index）/2］范圍內檢測最值點，該最值點所在位置即為最終的 A 點所在位置（符號記為 A_index）。
（10）最終的 B 點在 c（A_index：max_index）中幅值等于 0.15·xC 位置處。
（11）檢測 X 點。在 d_d［（1 + n·sp）：（n + 1）· sp］上，C 點右側一定范圍內檢測 X 點。考慮到算法效率，縮小 X 點的檢測范圍為 d_d{［2·min_index – QC_index –（min_index – QC_index）/2］：［2·min_index – QC_index +（min_index – QC_index）/2）］}，在該范圍內檢測距離 min_index 最近的一個過零點，對應位置為準 X 點所在位置（符號記為：QX_index）。
（12）判斷準 X 點 c［QX_index，x（QX_index）］是否為極值點，若是，則此點即為 X 點；否則，在 c［QX_index –（QX_index – min_index）/2：QX_index + （QX_index – min_index）/2］范圍內檢測極值點，將該極值點作為最終的 X 點。
（13）n = n + 1。重復步驟 （1）～（11），直到檢完整段信號。
2 實驗及結果
2.1 實驗數據來源
本次實驗中所用數據由重慶市武警重慶總隊醫院的心功能科室提供，數據的采集設備為心功能無創檢測分析儀（KF_ICG-303 型，重慶科發醫療器械有限公司，重慶）[21]。本文選取 30 例形態差異較明顯的患者的信號，對本文所提算法的有效性進行評估。其中，30 例患者的年齡分布于 31～88 歲，男性患者 17 例，女性患者 13 例。針對各例患者，分別選取時長約 5 min 的信號，對其特征信息進行識別、統計。
2.2 實驗結果及分析
本文隨機選取一例患者信號（采樣頻率為 1 000 Hz），對實驗過程及結果進行詳細論述如下。
部分原始信號如圖 2 所示。
圖2
部分原始信號
Figure2.
Unprocessed signal fraction
該信號夾雜噪聲，會影響對信號的特征提取結果，故首先對信號進行降噪處理。對信號進行 AEEMD 分解，AEEMD 分解后的各階 IMF 如圖 3 所示。
圖3
待處理信號 AEEMD 分解結果
Figure3.
Decomposed results in AEEMD of unprocessed signal
選取較高頻噪聲占主要成分的 IMFs4 信號進行疊加，并用 WT 法對其進行降噪處理，對基線漂移占主要成分的 IMF9、IMF10、IMF11 直接舍去，重構 WT 算法降噪后的 IMFs4 信號及有用信號占主要成分的 IMF5、IMF6、IMF7、IMF8 信號，最終對 ICG 信號的降噪結果如圖 4 所示。
圖4
降噪后信號
Figure4.
The de-noised signal
濾波前后，信號的頻譜變化如圖 5 所示。
圖5
信號降噪前后信號頻譜變化
Figure5.
The frequency spectrum of signal before and after de-noise
整體來看，降噪后信號中較高頻噪聲均被有效濾除，根據降噪前后信號的 20 Hz 以內的頻譜分布情況，較低頻噪聲和部分重疊頻率噪聲對應位置的頻率也得以明顯衰減，證明了 AEEMD_WT 法對 ICG 信號中多種類型的噪聲濾除的有效性。
此外，經反復實驗論證，在 IMF6、IMF7 疊加完進行差分后的信號中，各節拍信號的過零點位置即為對應 ICG 信號中的 C 點位置，故將此作為特征波檢測層。特征檢測層和降噪后信號的特征對應關系如圖 6 所示。
圖6
檢測層和降噪后信號特征對應關系
Figure6.
The corresponding relationship between the de-noised signal and the detection signal
如圖 6 所示，檢測層差分信號 d_d（n）中的橫虛線代表該信號的過零點位置，豎虛線表示了 d_d（n）和降噪后的 ICG 信號 c（n）的特征點的對應關系。在單個分段長度內，c（n）中 C 點位于 d_d（n）中的最大值（極大值）和最小值（極小值）間過零點的位置，本文選擇在 c（n）上的該極值對范圍內檢測極值點，以此確定 C_index；在 d_d（n）信號中，C_index 左側第一個過零點位置應為理想情況下的 A_index，C_index 右側第一個過零點位置為理想情況下 X_index，為避免噪聲的干擾，在 c（n）中，以 d_d（n）信號中的極大值、極小值位置為參考，確立 A 點和 X 點的檢測范圍，檢測對應的極值點以避免誤檢。特征點檢測效果如圖 7 所示。
圖7
特征識別結果
Figure7.
The result of characteristic recognition
參照標準 ICG 信號各特征點的位置關系，并經過和有多年相關臨床經驗的武警重慶總隊醫院的心功能科室醫生反復討論分析，上述 ICG 信號的各特征點均可認為被有效檢出。研究結果顯示，采用本文所提算法對臨床上大量不同病癥類型的信號進行分析處理，均取得了滿意的結果。最終的特征點檢測統計結果如表 1 所示。需要指出的是，本文所提出的特征點識別算法中，A 點、B 點、X 點的檢測均有賴于 C 點的成功檢出，即任一節拍中 C 點出現漏檢或誤檢，則整個節拍都將漏檢或誤檢；所謂誤檢，為錯檢和多檢的節拍總數。考慮到各特征信息的位置關系及基于此的算法設計原理，統計時，本文選擇以整個節拍為單位（任一特征點出現誤檢，即視為此節拍誤檢），參照公式（5）、（6）、（7），分別計算漏檢率、誤檢率和準確率，最終的統計結果如表 1 所示。
表1
特征點檢測結果統計
Table1.
The statistical result of feature points detected
|
|
|
特征檢測過程中，多處漏檢和誤檢均發生在信號自身特征缺失處。分析原因可能是由于 ICG 信號的微弱性使其在采集過程中極易受到噪聲的干擾，特別是噪聲頻率重疊于有用信號頻率的情況，導致信號特征被淹沒；另一方面，受試者在測量過程中偶見的保護性反射行為，如咳嗽等，均會影響信號形態特征；此外，受試者均有不同類型的心血管疾病，疾病的嚴重程度和類型會影響 ICG 信號的形態特征，如幅值驟減、信號變形、甚至部分特征點消失等。
3 結論
臨床上，基于胸部生物阻抗技術的心臟血流動力學分析的可靠性嚴重受限于 ICG 信號特征的準確識別，為保證心功能無創檢測儀更好地發揮其臨床應用價值，本文首先考慮到 ICG 信號所含噪聲的多樣性，提出了一種 AEEMD_WT 降噪算法，成功濾除了信號中較高頻噪聲、低頻噪聲及肌電干擾、電極移動等造成的和有用信號頻率有重疊的噪聲，為特征點的準確識別奠定了基礎；繼而，考慮到病理性 ICG 信號的形態多變性，設計了可自適應更新檢測節拍長度的 AEEMD_D_AS 特征檢測算法完成了對 A 點、B 點、C 點和 X 點的檢測；最后，以臨床上隨機抽樣到的 30 例病理性 ICG 信號對降噪及特征識別算法進行驗證，最終檢測的準確率高達 99.72%，進一步保證了基于胸部生物阻抗技術的心臟血流動力學參數的計算準確度。
