岳大超，甘良志，劉海寬，余南南

江蘇師范大學(xué) 電氣工程及自動化學(xué)院，江蘇 徐州 221116

1 引言

長久以來，心血管疾病一直是威脅人類健康的主要問題。目前，心血管病死亡占城鄉(xiāng)居民總死亡原因的首位，預(yù)計到2030年大約要有2 330萬人死于心血管病[1]。面對這種趨勢，對心血管病的早期診斷與預(yù)防顯得尤為重要。

傳統(tǒng)的診斷方法是患者到醫(yī)院，醫(yī)師對其進行心電圖檢查，進而給出診斷結(jié)果[2]，任務(wù)繁重且需要醫(yī)師有豐富的臨床經(jīng)驗和專業(yè)知識。而稀缺的醫(yī)療資源，很難滿足數(shù)量龐大患者群體的要求。為了提高就診效率，方便性和快捷性，出現(xiàn)了自動化診斷技術(shù)，輔助醫(yī)師進行診斷。

心率變異是指心動間期之間的時間變異數(shù)，其研究對象是心動間期。人的心率不是一成不變的，兩次心搏之間存在著微小的時間差異，計算心動間期的差異，即可了解心率變異性（Heart Rate Variability，HRV）。心率變異性可以評估心臟交感神經(jīng)與副交感神經(jīng)對心血管活動的影響，蘊含著心血管方面的大量信息[3-4]。臨床研究表明，心率變異性的降低是心肌梗死、高血壓、心絞痛等心血管疾病發(fā)病的癥狀。因此，心率變異性的研究，在評價心血管系統(tǒng)功能、預(yù)測心血管疾病的發(fā)作，以及為心血管疾病的早期診斷具有重要的意義[5-6]。

Poincare散點圖是心率變異性一種重要的研究方法。通過使用連續(xù)的心搏間期在直角坐標(biāo)系中繪制圖形，反映相鄰心搏間期的變化，能顯示心搏間期的特征。Poincare散點圖有多種形態(tài)，包括彗星狀、扇形等，不同的形狀反映不同的心臟狀態(tài)[7-8]。雖然Poincare散點圖是一種有效的心率變異性分析方法，但是并不能體現(xiàn)其隨時間變化的趨勢，對于某些心血管疾病、身體健康狀況等不能很好地體現(xiàn)其心率變異性性質(zhì)。于是，一些學(xué)者提出了改進策略，即一階差分散點圖，通過相鄰心搏間期的差值來繪制散點圖。然而，這種方法又丟失了原有的心搏間期絕對值信息。因此，在文獻[5]中，將二者結(jié)合起來，提出了一種RdR散點圖，以此來同時反映心搏間期及其變化。

目前，對于心率變異性分析的研究很多[9-12]。但是，并沒有根據(jù)散點圖來有效識別、區(qū)分不同的心血管疾病，實現(xiàn)心電自動化診斷的研究。本文即是通過簡化粒子群算法（Simplified Particle Swarm Optimization，SPSO）與集成稀疏核主分量分析（Integration Sparse Kernel Principal Component Analysis，ISKPCA）方法來對RdR散點圖進行識別分類，為心電自動化診斷的實現(xiàn)提供一些基礎(chǔ)。

2 簡化粒子群算法

基本粒子群算法（Particle Swarm Optimization，PSO）是一種智能優(yōu)化算法，算法首先在搜索范圍內(nèi)隨機初始化一群粒子，之后不斷迭代，迭代過程中更新兩個參數(shù)：一是粒子自身的最優(yōu)位置；一是種群的最優(yōu)位置；直至達到終止條件[13]。

胡旺等人把基本粒子群算法進行了簡化，稱為簡化粒子群算法（Simplified Particle Swarm Optimization，SPSO）。假設(shè)在D維搜索空間中，設(shè)置S個粒子，迭代T次，第i個粒子在d維的位置用Xid表示，用 pBestid記錄該粒子最好的位置，用gBestd記錄種群最好的位置。簡化后公式為：

式（1）中，右邊第一項是“粒子慣性”，為個體歷史對現(xiàn)在的影響，利用慣性權(quán)重ω來調(diào)節(jié)影響；第二項是“自我認(rèn)知”，是粒子對自身的思考，c1是學(xué)習(xí)因子；第三項是“社會引導(dǎo)”，表示粒子對歷史上最優(yōu)粒子的模仿，c2是學(xué)習(xí)因子，r1，r2是服從U(0,1)分布的隨機數(shù)。該算法僅通過粒子的位置進行更新迭代，把方程從二階降為了一階，并通過理論和實驗證明了該方法的有效性、優(yōu)越性，這種改進使得整個過程變得更為簡單。

3 集成稀疏核主分量分析

主分量分析是一種典型的無監(jiān)督算法，常用于解決原始空間的線性問題[14-15]，而為了在特征空間中用線性方法解決原始空間的非線性問題，B.Scholkopf等人提出了核主分量分析（Kernel Principal Component Analysis，KPCA）[16]。

其中，矩陣K是N×N的矩陣，也稱核矩陣。則式（3）問題轉(zhuǎn)換為：

過程可以直接在K上運算：

本文選擇一種常用的核函數(shù)，徑向基核函數(shù)來進行計算。當(dāng)選擇徑向基核函數(shù)時，會有明顯的過學(xué)習(xí)問題，主要是由于徑向基核函數(shù)對應(yīng)的特征空間是無限維的，通過該方法得出的主分量的次數(shù)與給定樣本的維數(shù)無關(guān)，而與樣本數(shù)量有關(guān)。為解決這個問題，引入了稀疏化方法，即為稀疏核主分量分析（Sparse Kernel Principal Component Analysis，SKPCA）。

在核主分量分析中，特征向量能用樣本{φ(x1),φ(x2),…,表示為 ，也就是說特征向量?是全部樣本的線性組合，這為特征向量的稀疏化提供了一種思路。通過近似基求解算法，來求{φ(x1),φ(x2),…,φ(xN)}的近似基，從而得到的稀疏化方法。

近似基求解的具體步驟是：

（1）建立集合Xl={φ(x1)}為近似最大無關(guān)組，XA=?。

（3）如果得到極小值value≤ε，則把對應(yīng)的元素加到XA當(dāng)中，否則為Xl。ε是線性相關(guān)截尾誤差，在有限的樣本中求無限維空間的基幾乎沒有稀疏性，因此選擇近似計算。

（4）返回步驟（2），直至完成所有的計算，計算結(jié)束。

其中，步驟（2）中的求解過程如下，令

可知式（11）是特征值與特征向量的問題，即

推導(dǎo)得：

其中，K(m,:)表示核矩陣K的第m行，K(:,n)表示第n列表示1行N 列的行向量。式（12）也就變成了：

為典型的特征值與特征向量問題。

雖然稀疏化可以控制學(xué)習(xí)機的學(xué)習(xí)能力，防止過學(xué)習(xí)，但是稀疏化方法仍然可能會丟失數(shù)據(jù)集的某些重要特性，即欠學(xué)習(xí)問題。因此，有必要引入集成方法，稱為集成稀疏核主分量分析（Integration Sparse Kernel Principal Component Analysis，ISKPCA）。由于核主分量分析是一種無監(jiān)督學(xué)習(xí)，不能從外部獲得指導(dǎo)，即不能對某次學(xué)習(xí)結(jié)果進行獎賞或懲罰，所以，本文選擇簡單平均的方法。SKPCA的集成方法具體步驟流程如下：

（1）設(shè)置重復(fù)的次數(shù)Re。（2）計算核矩陣K=ΦTΦ。

4 Rd R散點圖

令△RRn=RRn-RRn-1，RRn為第n個心搏間期。RRn為橫坐標(biāo)，△RRn為縱坐標(biāo)，逐次繪點，即得RdR散點圖。為驗證本文算法的有效性，使用的數(shù)據(jù)為MIT-BIH數(shù)據(jù)庫中的心電數(shù)據(jù)作為測試樣本，實驗過程中用到的數(shù)據(jù)編號如表1所示。

實驗當(dāng)中，使用MATLAB繪制來繪制RdR散點圖，取編號100中的部分?jǐn)?shù)據(jù)來繪制RdR散點圖，示例如圖1所示。

表1 本文使用的MIT-BIH數(shù)據(jù)庫的數(shù)據(jù)

圖1 編號100的RdR散點圖

5 實驗

將ISKPCA用于RdR散點圖的分類識別，需要先建立ISKPCA模型，然后選擇控制指標(biāo)。平方預(yù)測誤差（Squared Prediction Error，SPE）是常用的控制指標(biāo)。

對于第i個樣本Xi，假設(shè)經(jīng)過SKPCA求得其前n個主分量為t1,t2,…,tn，對應(yīng)的特征值為 λ1,λ2,…,λn。用表示特征空間N個主成分的重構(gòu)向量同理則樣本 X 的SPE定義為：

接下來要獲取樣本，將表1中所示的23組數(shù)據(jù)使用MATLAB繪制成RdR散點圖。在每個散點圖中，選擇1 000個RR間期來繪制散點圖，即至少10分鐘的數(shù)據(jù)，依據(jù)MIT-BIH的360 Hz采樣頻率計算，每幅散點圖大約需要使用216 000個心電數(shù)據(jù)點。為了盡可能多地獲得訓(xùn)練樣本，采用每次移動10個RR間期來繪制散點圖，共獲得樣本2 819例。算法具體步驟如下：

（1）讀取圖像數(shù)據(jù)，對圖像進行縮放到同一規(guī)格，轉(zhuǎn)成灰度圖，將一幅圖的多維矩陣數(shù)據(jù)轉(zhuǎn)換成一維，并對樣本進行分類標(biāo)記，分別用1、2、3、4、5對應(yīng)著正常竇性心律、非偶聯(lián)早搏、二聯(lián)律早搏、三聯(lián)律早搏以及混合早搏。

（2）對數(shù)據(jù)進行Z-score標(biāo)準(zhǔn)化處理。公式如下：

其中，μ、σ分別是數(shù)據(jù)的均值與方差。

（3）分別對每類樣本采樣，隨機抽取其中85%的數(shù)據(jù)作為訓(xùn)練樣本，剩余的15%作為測試樣本。

（4）使用簡化粒子群算法搜尋參數(shù)，包括求解近似基的誤差參數(shù)ε、高斯核函數(shù)參數(shù)gamma以及核主分量分析控制限t的值，以準(zhǔn)確率為適應(yīng)度函數(shù)，分別對每類樣本求解近似基、特征值、特征向量。

（5）分別計算每個測試樣本的SPE，其值與某類樣本的SPE差值最小者為測試樣本的預(yù)測類別，與實際類別比較，計算準(zhǔn)確率，若滿足要求則步驟（6），反之返回步驟（4）。

（6）使用所得分類模型，進行樣本交叉實驗，測試算法性能。

本文分類模型參數(shù)最終選擇為：ε=250，高斯核函數(shù)參數(shù)gamma=600，主分量控制t=90%，集成個數(shù)Re=10。經(jīng)過多次樣本交叉驗證，使用KPCA平均識別率為72%，SKPCA的平均識別率約為91%，集成SKPCA的平均識別率為98.6%。與文獻[17]中使用BP神經(jīng)網(wǎng)絡(luò)來實現(xiàn)3種心律識別的97%的正確率相比，本文算法準(zhǔn)確率有所提高。由此可以看出，RdR散點圖可以反映上述幾種心律的特征，另外，也體現(xiàn)了使用本文算法來進行RdR散點圖分類識別的有效性。

本文使用簡化粒子群算法的參數(shù)：種群規(guī)模是20，進化100代，搜索空間依據(jù)實驗測試設(shè)置在[0，1 000]。從測試樣本中隨機選擇某一類的心電數(shù)據(jù)樣本，如圖2選擇到實際類別為“3”的一些測試樣本集，將其分別與各個類別計算其SPE的差值，結(jié)果如圖2，可以看到其與類別為“3”的SPE差值最小，分類正確；圖3是隨機抽取部分樣本的預(yù)測結(jié)果與實際標(biāo)記樣本對比，可以看到全部分類正確。

圖2 測試樣本與各類的SPE差值

6 結(jié)束語

圖3 部分預(yù)測結(jié)果與實際樣本對比

核主分量分析存在著過學(xué)習(xí)的問題，根據(jù)統(tǒng)計學(xué)習(xí)理論可知，稀疏化方法可以控制學(xué)習(xí)機的學(xué)習(xí)能力，防止過學(xué)習(xí)，而且通過稀疏化可以削減樣本數(shù)量，減少計算量，而為了應(yīng)對由于稀疏化使某些特征被遺漏的問題，即欠學(xué)習(xí)，又引入了集成方法。本文用MIT-BIH中的心電數(shù)據(jù)來驗證本文算法的有效性，算法使用心電數(shù)據(jù)繪制RdR散點圖，通過簡化粒子群算法對集成稀疏核主分量分析的參數(shù)進行尋優(yōu)，以正確識別率為適應(yīng)度函數(shù)，并用所學(xué)習(xí)到的分類模型來對五種心律的RdR散點圖進行分類識別，實驗結(jié)果顯示，分類效果較好。但一方面，由于本文使用的數(shù)據(jù)樣本及類別相對于龐大的心血管疾病種類而言都較少，使得不能非常全面、客觀地反映所有問題；另一方面，本文所述方法，步驟較為繁瑣，計算周期長。接下來的研究重點是：收集更多的心電樣本數(shù)據(jù)；尋找更簡單，更有效的方法、更高效的方法。為了實現(xiàn)心電自動化診斷還需要進一步的研究。