心率變異性多參數(shù)分析的BP網(wǎng)絡(luò)用于心衰診斷的研究

侯鳳貞 黃曉林 寧新寶*

1(南京大學(xué)電子科學(xué)與工程學(xué)院生物醫(yī)學(xué)電子研究所，南京 210093)

2(中國藥科大學(xué)基礎(chǔ)部信息管理教研室，南京 210009)

引言

心臟病是一種常見的多發(fā)慢性疾病，近年來，隨著生活環(huán)境的變更、生存節(jié)奏的加快以及社會(huì)競爭和各方面壓力的加劇，心臟病的患病率逐年升高。在我國，每年有幾十萬人死于心臟病，全世界更有1/3人口的死亡是由心臟病引起的。因此，如何有效地檢測與評(píng)價(jià)心臟的功能狀況，以便對(duì)心臟疾病進(jìn)行準(zhǔn)確的預(yù)報(bào)和診斷，是目前治療心臟疾病的一個(gè)重要研究課題［1］。

在臨床診斷心臟疾病時(shí)，作為一種無創(chuàng)檢測手段，體表心跳頻率變異信號(hào)(Heart Rate Variability，HRV)分析受到很大的重視。目前對(duì)于 HRV信號(hào)的分析，常采用時(shí)域分析法、頻域分析法以及非線性動(dòng)力學(xué)分析法［1－2］。

時(shí)域分析方法是研究HRV的最簡單易實(shí)現(xiàn)的方法，主要是借用統(tǒng)計(jì)學(xué)的方法，對(duì)時(shí)間序列進(jìn)行統(tǒng)計(jì)性分析。研究表明，在冠心病患者，嚴(yán)重的心肌病患者、充血性心臟病、心肌梗塞患者及老年患者中整個(gè)HRV時(shí)域分析指標(biāo)的下降預(yù)示著術(shù)后狀態(tài)較差和死亡威脅的增加［1］。時(shí)域分析已經(jīng)廣泛的應(yīng)用來比較心肌梗死的治療效果，估計(jì)冠狀動(dòng)脈受阻，預(yù)測心律失常等等［1］。

頻域分析法即功率譜分析法，是將一定時(shí)間內(nèi)的連續(xù)RR間期序列進(jìn)行離散傅里葉變換得到以頻率為橫坐標(biāo)、功率為縱坐標(biāo)的功率譜曲線的方法。頻域分析法提供了能量隨頻率變化分布的基本信息，能夠較為有效地鑒別出交感神經(jīng)和迷走神經(jīng)各自對(duì) HRV 的影響［2］。

近些年來，人們?cè)絹碓角宄卣J(rèn)識(shí)到，心臟電活動(dòng)與傳統(tǒng)研究的隨機(jī)信號(hào)和周期信號(hào)都有所不同，它由多種獨(dú)立因素協(xié)調(diào)控制，猶如一個(gè)遠(yuǎn)離平衡點(diǎn)的多輸入物理系統(tǒng)，富含復(fù)雜的非線性成分［1］。因此，人們開始廣泛應(yīng)用非線性的概念和方法對(duì)心臟動(dòng)力學(xué)進(jìn)行研究。眾多的非線性分析方法，如冪律分析、熵分析等已經(jīng)廣泛地用于HRV信號(hào)的分析中;許多衡量HRV非線性特性的指標(biāo)，如尺度指數(shù)、樣本熵值、李雅普諾夫指數(shù)等被陸續(xù)提了出來［2－5］。另外，在對(duì) HRV信號(hào)的分析中結(jié)合一些行之有效的研究手段，如多尺度分析、符號(hào)動(dòng)力學(xué)分析、時(shí)間逆轉(zhuǎn)分析，也都獲得了一定的進(jìn)展［6－8］。

但是，對(duì)照比較充盈性心衰(congestive heart failure，CHF)患者和正常人的各指標(biāo)，我們發(fā)現(xiàn)很難用一個(gè)單一的指標(biāo)來完全區(qū)分CHF患者和健康人。鑒于單一指標(biāo)所表達(dá)出來的信息具有片面性，且各指標(biāo)的使用都有其優(yōu)缺點(diǎn)，聯(lián)合多個(gè)指標(biāo)進(jìn)行分析有可能提高CHF的診斷正確率。作為生物控制論的一個(gè)顯著成果，人工神經(jīng)網(wǎng)絡(luò)的工作原理和功能特點(diǎn)接近于人腦，不是按照給定的程序一步步地機(jī)械執(zhí)行，而是根據(jù)自身適應(yīng)環(huán)境、總結(jié)規(guī)律、完成運(yùn)算、識(shí)別和控制工作。目前，人工神經(jīng)網(wǎng)絡(luò)在模式識(shí)別與分類領(lǐng)域中已經(jīng)取得了較為顯著的成果［9］。因此在本研究選取了多個(gè)從 HRV信號(hào)中提取的指標(biāo)，聯(lián)合起來作為人工神經(jīng)網(wǎng)絡(luò)的特征參數(shù)，利用誤差反向傳播(back propagation，BP)算法，建立起用于心衰診斷的智能模型，擬為充盈性心衰的診斷提供了一條無創(chuàng)、簡捷、高準(zhǔn)確率的補(bǔ)充途徑。

1 實(shí)驗(yàn)數(shù)據(jù)與方法

1.1 實(shí)驗(yàn)數(shù)據(jù)

復(fù)雜生理信號(hào)研究資源庫 PhysioNet是由NIBIB(National Institute of Biomedical Imaging and Bioengineering)和 NIGMS(NationalInstituteof General Medical Sciences)聯(lián)合贊助的開放服務(wù)資源庫［10］。本研究使用的數(shù)據(jù)來源于其中的充盈性心衰chfdb和chf2db數(shù)據(jù)庫，以及正常竇性心率數(shù)據(jù)庫nsrdb和nsr2db。其中包括72例健康人，平均年齡:(54.6±16.0)歲，44例 CHF患者，平均年齡(55.6±11.3)歲。幾個(gè)數(shù)據(jù)庫中樣本采樣時(shí)間都在18 ～24 h之間，為長時(shí)采樣。但是在計(jì)算各HRV指標(biāo)時(shí)，考慮到晝夜節(jié)律的影響，數(shù)據(jù)長度都選定為20000點(diǎn)(約5 h的采樣時(shí)間)。

1.2 常用HRV指標(biāo)的計(jì)算

1.2.1 時(shí)域指標(biāo)的計(jì)算

對(duì)于數(shù)據(jù)長度為N的HRV序列{RRi:1≤i≤N}，借用統(tǒng)計(jì)學(xué)的方法進(jìn)行統(tǒng)計(jì)性分析性時(shí)可以得到時(shí)域分析中的長時(shí)指標(biāo)，常用的指標(biāo)及其計(jì)算方法如下［1］:

· 全部正常竇性RR間期均值aRR，可依下式計(jì)算得到:

·全部正常竇性RR間期標(biāo)準(zhǔn)差SDNN，可依下式計(jì)算得到:

·全程相鄰RR間期之差的均方根RMSSD，可依下式計(jì)算得到:

·每5 minRR間期平均值的標(biāo)準(zhǔn)差SDANN和每5 min的RR間期標(biāo)準(zhǔn)差的平均值SDNN index:本研究中將每個(gè)長度為N的HRV序列按每段500點(diǎn)劃分成K=N/500小段，先計(jì)算每小段內(nèi)的RR間期平均值，再計(jì)算這K個(gè)平均值的標(biāo)準(zhǔn)差得到 SDANN指標(biāo);而若先計(jì)算每小段內(nèi)RR間期的標(biāo)準(zhǔn)差，再計(jì)算這K個(gè)標(biāo)準(zhǔn)差的平均值則得到SDNNindex指標(biāo)。

1.2.2 頻域指標(biāo)的計(jì)算

按照國際上的標(biāo)準(zhǔn)，將HRV序列進(jìn)行離散傅里葉變換，并對(duì)得到的功率譜依據(jù)頻率的不同可以劃分為4個(gè)頻段［2］:0～0.003 Hz的超低頻段(ultra low frequency，ULF)、0.003～0.04 Hz的極低頻段(very low frequency，VLF)、0.04～0.15 Hz的低頻段(low frequency，LF)和 0.15～0.4 Hz的高頻段(high frequency，HF)。

對(duì)于長時(shí)HRV信號(hào)的頻域分析，通常用的指標(biāo)有上述4個(gè)頻段的功率［1～2］，在本研究中我們分別記為 ULFP(ULF Power)、VLFP(VLF Power)、LFP(LF Power)、HFP(HF Power)。此外，所有小于0.4 Hz的頻率成分的總功率TP(Total Power)也是一個(gè)常用指標(biāo)［1］。

1.2.3 非線性指標(biāo)的計(jì)算

在本研究中，對(duì)數(shù)據(jù)長度為 N的 HRV序列{RRi:1≤i≤N}，計(jì)算如下幾個(gè)非線性指標(biāo)，這些指標(biāo)最近被報(bào)道能較成功地區(qū)分CHF患者和健康人:

·MSE4:是對(duì)RRi的增量序列做多尺度樣本熵分析時(shí)，在尺度4下的熵值［11］。其計(jì)算方法可依如下步驟進(jìn)行。

步驟1:獲得增量序列

步驟2:對(duì)上述增量序列進(jìn)行尺度為4的粗?；?

y序列的長度為K=(N－1)/4。

步驟3:對(duì)序列y做m維嵌入，得到矢量序列:B(m)(i)=(yi，yi+1，…，yi+m－1)，1 ≤ i≤ K － m+1

步驟4:當(dāng)兩個(gè)矢量 B(m)(i)和 B(m)(j)的距離小于一個(gè)設(shè)定值r時(shí)，則認(rèn)為兩矢量是相同的。定義n(m)i是與B(m)(i)相同的矢量個(gè)數(shù)，則m維嵌入下的矢量相同的概率為:

步驟5:將嵌入維數(shù)增加到 m+1，同法計(jì)算C(m+1)(r).

步驟6:則y序列的樣本熵值，也就是MSE4指標(biāo)值為

·P8:是對(duì) RRi做高維的時(shí)間不可逆性分析，在嵌入維為8時(shí)參數(shù)Pm的值［12］。其計(jì)算方法包括以下步驟。

步驟 1:對(duì) RRi序列做 8維嵌入，得到矢量序列:

步驟2:計(jì)算該矢量序列在各個(gè)二維平面(RRi，RRi+n)(1≤n≤7)上的投影關(guān)于該平面主對(duì)角線RRi=RRi+n的不對(duì)性測度:

式中，ΔRRn=RRi+n－ RRi，N(ΔRRn)表示值小于 0的ΔRRn的個(gè)數(shù)，N(ΔRRn≠0)則表示值不等于0的ΔRRn的個(gè)數(shù)。

步驟3:所求的非線性指標(biāo)值為

· POLVAR20:這是一個(gè)符號(hào)動(dòng)力學(xué)參數(shù)。其計(jì)算方法包括以下步驟［13］。

步驟4:對(duì) RRi的增量序列 ΔRRi(1≤i≤N－1)，按照絕對(duì)值小于20 ms符號(hào)化為0，否則為1的方式，符號(hào)化僅含0、1兩個(gè)符號(hào)的序列Di(1≤i≤N－1);

步驟5:對(duì)Di序列做6維嵌入，得到矢量序列:

步驟6:統(tǒng)計(jì)矢量“000000”在矢量序列 B(i)中出現(xiàn)的次數(shù)n，即可得到所求的指標(biāo)值:

1.3 HRV指標(biāo)的評(píng)價(jià)

一般以受試者工作特征曲線(receiver operator characteristic curve，ROC)下的面積(area under the curve，AUC)為依據(jù)來判斷指標(biāo)的診斷價(jià)值［14］。如果定義靈敏度(sentivity，sen)為真陽性率，即被正確診斷為異常的人數(shù)與總異常人數(shù)的比值，定義特意度(specificity，spe)為真陰性率，即被正確診斷為正常的人數(shù)與總正常人數(shù)的比值，那么以(1－spe)為橫軸，以sen為縱軸，將待檢驗(yàn)的指標(biāo)值按不同的閾值用最基本的線性感知分類器分類，分別計(jì)算 sen和spe，將結(jié)果繪制于圖上，就可以得到 ROC曲線。一般AUC值應(yīng)在1.0和0.5之間。在 AUC＞0.5的情況下，AUC越接近于 1，說明診斷效果越好。AUC在0.5～0.7時(shí)有較低準(zhǔn)確性，AUC在 0.7～0.9時(shí)有一定準(zhǔn)確性，AUC在 0.9以上時(shí)有較高準(zhǔn)確性。AUC=0.5時(shí)，說明診斷方法完全不起作用，無診斷價(jià)值［14］。

1.4 BP網(wǎng)絡(luò)的設(shè)計(jì)

三層的 BP網(wǎng)絡(luò)(即一個(gè)輸入層、一個(gè)隱層、一個(gè)輸出層)可以擬合任意的非線性函數(shù)［15］。在本研究中，為了聯(lián)合時(shí)域、頻域、非線性動(dòng)力學(xué)分析方法，選擇常用的時(shí)域指標(biāo)中AUC值最大的指標(biāo)，頻域指標(biāo)中AUC值最大的指標(biāo)以及3個(gè)非線性指標(biāo)MSE4、P8、POLVAR20，共同作為 BP 網(wǎng)絡(luò)的特征向量，因此輸入層神經(jīng)元即為5個(gè)。網(wǎng)絡(luò)的輸出層神經(jīng)元個(gè)數(shù)K由需要分類的種數(shù)決定，本研究只研究健康人和CHF患者的區(qū)別，故選定輸出神經(jīng)元為一個(gè)。而為了不增加網(wǎng)絡(luò)的復(fù)雜程度，確定網(wǎng)絡(luò)隱層為一層，隱層神經(jīng)元數(shù)定為7。隱層傳遞函數(shù)選擇非線性激活函數(shù)tansig();輸出層傳遞函數(shù)亦選擇非線性激活函數(shù)tansig()。

1.5 BP網(wǎng)絡(luò)的訓(xùn)練與評(píng)估

采用 Levenberg－Marquard 算法［16］來訓(xùn)練如上的BP網(wǎng)絡(luò);同時(shí)為了提高網(wǎng)絡(luò)的泛化能力，采用及早停止法(early stopping)來控制訓(xùn)練過程［17］。在這種控制模式下，網(wǎng)絡(luò)開始訓(xùn)練前需將全部的116個(gè)樣本按照6∶2∶2的比例隨機(jī)劃分為訓(xùn)練(Train)樣本、驗(yàn)證(Validation)樣本、測試(Test)樣本3個(gè)子集。使用訓(xùn)練樣本進(jìn)行學(xué)習(xí)，同時(shí)監(jiān)測驗(yàn)證樣本的誤差。驗(yàn)證樣本的誤差在網(wǎng)絡(luò)訓(xùn)練的初期是逐漸減小的，而當(dāng)網(wǎng)絡(luò)出現(xiàn)過擬合時(shí)，其誤差便開始上升，此時(shí)就應(yīng)停止訓(xùn)練。測試樣本用于比較不同的網(wǎng)絡(luò)模型，它的誤差開始增加的時(shí)間應(yīng)接近驗(yàn)證樣本。

考慮到一次訓(xùn)練的結(jié)果可能并不具有代表性，因此，本研究中共進(jìn)行了10000次的隨機(jī)劃分訓(xùn)練、驗(yàn)證、測試樣本的訓(xùn)練。在每一次的訓(xùn)練結(jié)束后，借助matlab神經(jīng)網(wǎng)絡(luò)工具箱中的sim函數(shù)，將全部的116個(gè)樣本的特征向量輸入至訓(xùn)練好的網(wǎng)絡(luò)進(jìn)行仿真。對(duì)于健康人樣本，期望網(wǎng)絡(luò)輸出值為-1，其仿真輸出值不大于0時(shí)視為判斷正確;對(duì)于CHF患者，期望輸出為1，其仿真輸出值大于0時(shí)視為判斷正確。由于每次訓(xùn)練時(shí)的測試樣本子集都是隨機(jī)確定的，可能每次的測試子集都不相同，因此，本研究中使用全樣本集的識(shí)別正確率，即116個(gè)樣本中判斷正確的樣本百分?jǐn)?shù)，而不是用測試子集的識(shí)別正確率來評(píng)估神經(jīng)網(wǎng)絡(luò)的建模結(jié)果。

2 結(jié)果

對(duì)116個(gè)HRV樣本數(shù)據(jù)分別進(jìn)行了時(shí)域、頻域、非線性動(dòng)力學(xué)分析，分別計(jì)算了上述各種常用指標(biāo)及各指標(biāo)在診斷CHF時(shí)的AUC值，結(jié)果如表1所示。從表1中選擇出AUC值最大的時(shí)域指標(biāo)SDNN、頻域指標(biāo) LFP、三個(gè)非線性指標(biāo) MSE4、P8、POLVAR20共同作為 BP網(wǎng)絡(luò)的輸入向量，訓(xùn)練10000次后，對(duì)于全樣本集的識(shí)別正確率平均為86.97%，標(biāo)準(zhǔn)差為8.86%，而最大識(shí)別正確率高達(dá)99.14%，即如圖1所示，116個(gè)樣本中僅有編號(hào)為100的樣本被誤判。而從圖2則可以看出僅經(jīng)過21步的迭代后，驗(yàn)證集的誤差開始呈現(xiàn)上升趨勢，網(wǎng)絡(luò)停止訓(xùn)練;測試集的誤差也在第22步迭代后開始增加，說明網(wǎng)絡(luò)已經(jīng)取得了較好的泛化能力。

表1 健康人組和CHF組常用HRV指標(biāo)值(均值 ±標(biāo)準(zhǔn)差)及 AUC值Tab.1 The values(mean±standard errors)of HRV indicesfor both healthy and CHF group and the corresponding AUC values

3 討論與結(jié)論

作為對(duì)比，選擇表1中AUC值最大的MSE4指標(biāo)，編程計(jì)算得出:僅僅采用該指標(biāo)進(jìn)行 CHF組和健康人組分類時(shí)，對(duì)于實(shí)驗(yàn)中的116個(gè)樣本的最大正確率為85.34%。

圖1 最大正確率時(shí)116個(gè)樣本的仿真值Fig.1 The simulation results for the 116 samples when the maximum accuracy achieved

圖2 最大正確率時(shí)網(wǎng)絡(luò)訓(xùn)練過程的平方誤差曲線Fig.2 The mean square error curve of network training when the maximum accuracy achieved

分別考察僅采用5個(gè)頻域參數(shù) ULFP、VLFP、LFP、HFP、TP，僅采用4個(gè)有一定診斷價(jià)值的時(shí)域參數(shù) aRR、SDNN、SDANN、SDNNindex，僅采用三個(gè)非線性參數(shù) MSE4、P8、POLVAR20，以及僅采用表 1 中AUC 值最大的 5個(gè)參數(shù) MSE4、VLFP、LFP、HF、TP 作為BP網(wǎng)絡(luò)的特征參數(shù)來對(duì)上述116個(gè)樣本進(jìn)行分類，同樣地每種試驗(yàn)都進(jìn)行了10000次的隨機(jī)劃分訓(xùn)練、驗(yàn)證、測試樣本的訓(xùn)練，并逐次進(jìn)行全部樣本的仿真測試，結(jié)果如表2所示。同時(shí)，t檢驗(yàn)的結(jié)果確證，聯(lián)合時(shí)域指標(biāo)SDNN、頻域指標(biāo)LTP、及3個(gè)非線性指標(biāo)時(shí)網(wǎng)絡(luò)的識(shí)別性能顯著高于僅上述采用某一類指標(biāo)時(shí)的識(shí)別性能(t檢驗(yàn)的P值均小于0.05)。

表2 采用不同特征參數(shù)的BP網(wǎng)絡(luò)用于CHF患者和健康人分類時(shí)的分類結(jié)果Tab.2 The classification results of BP neural network with different features used in CHF and healthy groups

上述對(duì)比實(shí)驗(yàn)的結(jié)果表明，聯(lián)合時(shí)域、頻域及非線性指標(biāo)作為神經(jīng)網(wǎng)絡(luò)的特征參數(shù)來區(qū)分兩類人群比僅使用某一個(gè)或者某一類指標(biāo)來分類，能獲得更高的識(shí)別正確率。因此可認(rèn)為，對(duì)于HRV信號(hào)的分析，無論是時(shí)域分析方法、頻域分析方法還是非線性動(dòng)力學(xué)分析方法，所得到的每一個(gè)指標(biāo)都只能片面地反映心臟狀態(tài)的一個(gè)方面，只有聯(lián)合線性(時(shí)域分析、頻域分析)和非線性方法，使用盡可能多的指標(biāo)才可能全面地揭示心臟的動(dòng)力學(xué)特征，從而提高心臟疾病的診斷正確率。雖然本研究中僅僅對(duì)心衰的診斷做了一定的探索，但從探索結(jié)果來看，聯(lián)合HRV信號(hào)分析的線性和非線性方法，使用多參數(shù)來進(jìn)行心臟疾病的診斷，極有可能取得較高的診斷正確率。因此，將進(jìn)一步充實(shí)臨床數(shù)據(jù)，繼續(xù)探討本方法在不同心臟疾病診斷中的應(yīng)用效果。同時(shí)，考慮到臨床上采集長時(shí)心電數(shù)據(jù)的不方便性，基于短時(shí)HRV信號(hào)分析的智能化診斷系統(tǒng)也是未來的一個(gè)研究重點(diǎn)。

［1］黃曉林.心率變異性的分析方法研究［D］.南京:南京大學(xué)，2008.

［2］Task Force.Heart rate variability:standard of measurement，physiological interpretation，and clinical use ［J］.European Heart Journal，1996，17(3):354 －381.

［3］Ivanov P C，Chen Zhi，Hu Kun，et al.Multiscale aspects of cardiac control［J］.Physica A，2004，344:685 － 704.

［4］Peng Chungkang，Havlin S，Stanley H E，et al.Quantification of scaling exponents and crossover phenomena in nonstationary heartbeat time series［J］.Chaos，1995，5(1):82 － 87.

［5］Kobayashi M，Musha T.1/f fluctuation of heartbeat period［J］.IEEE Transactions on Biomedical Engineering，1982，29:456 －457.

［6］Costa M，Goldberger AL，Peng Chungkang.Multiscale entropy analysis of complex physiologic time series［J］.Physical Review Letters，2002，89:068102.

［7］Li Jin，Ning Xinbao.Dynamic complexity detection in shortterm physiological series using base－scale entropy［J］.Physical Review E，2006，73:052902.

［8］Hou Fengzhen，Zhuang Jianjun，Bian Chunhua，et al.Analysis of heartbeat asymmetry based on multi－scale time irreversibility test［J］.Physica A，2010，389:754 －760.

［9］王瑞，李杰平，盧志剛，等.基于人工神經(jīng)網(wǎng)絡(luò)遙感圖像分類的應(yīng)用研究［J］.科技情報(bào)開發(fā)與經(jīng)濟(jì)，2011，21(3):108－111.

［10］Goldberger AL，Amaral LAN，Glass L，et al.PhysioBank，PhysioToolkit，and PhysioNet:components of a new research resource for complex physiologic signals ［J］. Circulation，2000，101(23):215－220.

［11］Huang Xiaolin，Ning Xinbao，Wang Xinlong.Multiscale entropy of heart beat interval increment series and its application in diagnosis of congestive heart failure ［J］. Chinese Science Bulletin，2009，54(20):3784－3789.

［12］Hou Fengzhen，Ning Xinbao，Zhuang Jianjun，et al.Highdimensionaltime irreversibility analysisofhuman interbeat intervals［J］.Medical Engineering and Physics，2011，33(5):633－637.

［13］Wessel N，Riedl M，Kurths J. Is the normalheartrate“chaotic”due to respiration?［J］CHAOS，2009，19:028508.

［14］顏虹.醫(yī)學(xué)統(tǒng)計(jì)學(xué)［M］.北京:人民衛(wèi)生出版社，2005.217－226.

［15］寧新寶，黃曉林，徐寅林，等.神經(jīng)網(wǎng)絡(luò)分析方法用于心臟病診斷的研究［J］.生物醫(yī)學(xué)工程學(xué)雜志，2004，21(2):284－287.

［16］趙弘，周瑞祥，林廷圻.基于 Levenberg－Marquardt算法的神經(jīng)網(wǎng)絡(luò)監(jiān)督控制［J］.西安交通大學(xué)學(xué)報(bào)，2002，36(5):523－527.

［17］戴榮，李仲奎.三維地應(yīng)力場BP反分析的改進(jìn)［J］.巖石力學(xué)與工程學(xué)報(bào)，2005，24(1):83－88.