張詩雨，楊 珂，夏春明＊＊，金陳玲，王憶勤，燕海霞

（1. 華東理工大學機械與動力工程學院 上海 200237；2. 上海中醫(yī)藥大學四診信息綜合實驗室 上海 201203）

1 引言

脈診作為傳統(tǒng)中醫(yī)四診的重要組成部分，其客觀化研究在中醫(yī)現(xiàn)代化發(fā)展中備受關注，基于不同原理的脈象采集裝置以及脈象信號分析處理也得到了較大發(fā)展。在脈象的分類和模式識別方面，對于傳統(tǒng)中醫(yī)上基本脈象類別的分類和某些疾病對應病脈的分類都已有較多研究，且根據(jù)時域分析、頻域分析、時-頻聯(lián)合分析等多種方法提取出大量特征用于脈象信號分類分析[1]。

Zhang DY 等人[2]對多普勒超聲裝置采集的橈動脈血流速度信號進行分析處理，根據(jù)脈象波形特征點提取了時域特征，并結合通過希爾伯特-黃變換(HHT)提取的頻譜特征，使用支持向量機(SVM)區(qū)分健康人群和膽囊炎以及腎炎患者，三類平均準確率達到75.9%。史紅斐[3]等人對脈象信號進行了五個尺度上的二進離散小波變換，將前四個頻帶的能量作為特征，使用SVM 對滑、弦、細、澀四類脈象進行分類，平均準確率達到87.5%。此外，還有研究人員提取了高斯混合模型參數(shù)特征[4，5]、小波包特征[6]、倒譜特征[7]等用于脈象的分類識別。這些特征結合起來包含了大量冗余性及差異性，Lei Liu 等人[8]對七種不同的特征提取方法（基于特征點的時域特征、自回歸模型特征、希爾伯特-黃變換特征、近似熵、小波包和小波變換）提取出的七類特征進行集成，使用多核學習SimpleMKL 算法對健康人群、糖尿病患者、腎病患者和胃病患者進行區(qū)分，最終平均識別準確率達到94%。雖然此研究對多種脈象特征進行了集成應用，但并未對特征進行分析評估。

目前大多應用的分類算法沒有對特征組合進行篩選分析，因此也難以充分利用脈象特征的冗余性和差異性，采用的特征數(shù)量過多不僅會增大計算量，而且會造成過擬合，降低分類正確率。為嘗試解決以上問題，本文對采集到的脈象信號進行濾波、歸一化等預處理后，提取脈象信號的時域、頻域及時頻特征共93維，使用集成學習中的隨機森林算法對特征的重要性進行排序，并使用SVM、BP-NN 以及隨機森林算法驗證重要性排序的正確性，確定應用于四類脈象分類的最佳脈象特征區(qū)間，據(jù)此對脈象特征進行降維，確定四類脈象適合的特征類別。

2 脈象數(shù)據(jù)預處理及特征提取

2.1 脈象信號預處理

脈搏波是由心臟持續(xù)不斷地跳動引起血液在血管中流動而造成的脈的搏動，當機體發(fā)生病變或臟腑失調時，血管、脈氣、脈血發(fā)生改變形成特定的“象”[9]。從頻域來看，脈象信號是一種周期性較強的準周期信號，也是時頻隨機變化的非平穩(wěn)信號，其頻譜成分主要分布于0-20 Hz之間[10]。由于脈象信號在采集過程中易受外界環(huán)境干擾，因此需要對采集到的數(shù)據(jù)進行去噪處理，并根據(jù)后續(xù)信號分析需求進行周期分割和特征提取。

2.1.1 脈象信號去噪

如表1所示，本文研究使用的脈象數(shù)據(jù)由上海中醫(yī)藥大學提供，包括滑脈、平脈、弦脈及實脈四種脈象共175例，診斷結果均由兩名中醫(yī)師同時做出并結論一致，采樣頻率為720 Hz，在最佳取脈壓力下記錄60 s得到。

本文采用小波變換進行脈象信號降噪，首先將含有噪聲的脈象信號進行多尺度小波分解，根據(jù)得到的高頻細節(jié)分量和低頻近似分量進行閾值量化處理，再經小波重構即可得到去噪后信號。Symlets 作為雙正交小波，具有有限緊支撐和近似對稱的特性，在尺度變換上較靈活且可降低信號重構的相移。使用小波基Sym8 對脈象信號進行10 層小波分解，在采樣頻率720 Hz 下，根據(jù)Nyquist 采樣定理，信號的頻率范圍為0-360 Hz，低頻分量 A10 的頻寬為 0-0.35 Hz，可用于去除由于人體移動而引起的基線漂移。高頻細節(jié)分量D1-D3頻寬為360-45 Hz包含了高頻噪聲和工頻干擾可采用軟閾值去噪處理。將經過閾值去噪后的1-3尺度上的細節(jié)分量和未處理的4-10 尺度細節(jié)分量進行小波重構，即可得到濾除高頻、工頻干擾和基線漂移的脈象信號，如圖1即為小波去噪后結果。

表1 本文用于分析的脈象信號類型及數(shù)量

2.1.2 脈象信號分割及平均

在提取脈象信號時域特征時，需要對一個脈搏周期內的特征點進行定位提取，由于一個受試者在連續(xù)采集得到的多周期信號之間會存在一定的差異，因此需要對同一個受試者的多周期脈象信號進行分割和平均。脈象信號的分割前需要找出多周期信號的谷峰對，通過設置距離閾值和高限閾值確保得到準確的波形起始點和峰值點。信號分割后將各個單周期的起始點和結束點之間的連線旋轉至水平位置再進行標準化，將標準化之后的單周期信號在相關性最大處進行對齊后求平均波形。對于各個樣本得到的平均波形可能存在數(shù)據(jù)長度的差異，為了使提取的時域特征更具備可比性，對平均后的單周期脈象信號進行重采樣，統(tǒng)一各樣本周期長度，如圖2 為經過周期分割、旋轉、標準化、平均和重采樣之后的脈搏波形。

2.2 脈象信號特征提取

為了研究脈象信號特征對四類脈象信號分類的影響，本文提取了脈象分類常使用的脈象時域、頻域以及時-頻域特征。

圖1 小波去噪前后對比圖

圖2 重采樣后單周期脈搏波形

2.2.1 時域特征提取

對每個樣本經過重采樣和平均后的單周期脈象信號進行時域特征提取，主要反映脈搏波形狀特點。如圖3 所示為一個單周期脈搏波的時域基本信息點，本文提取了28 種時域特征[11]，包括5 個幅值比例特征以及主波高度 2/3處的橫坐標寬度w），5 個時間比例特征，10 個斜率相關特征（升支斜率，降支在主波與降中峽間的斜率，降中峽至主波終點的斜率，主波升支上升最快點的斜率值及坐標值，兩個降支（lbc和lde）中下降最快點的斜率值及坐標值，以及上升最快點與下降最快點之間的橫坐標距離），7 個面積相關特征（收縮期脈圖面積As(曲線段abcde下)，舒張期脈圖面積Ad(曲線段efg下)和脈圖面積特征量脈搏波起點到主波峰之間的脈搏波波形與這兩點間的直線所圍成的面積Asmp，脈搏波主峰到降中峽之間的波形與這兩點間直線所圍成的面積Ampv，降中峽到重搏波峰間脈搏波與這兩點的直線所組成的面積Avsp，以及重博波峰到波形末尾間脈搏波與這兩點的直線所組成的面積Aspe[12]），1 個脈搏波類型相關特征（如圖4 為四種脈搏波類型，其中type值即代表了樣本脈搏波類型）。

圖3 單周期脈搏波時域特征點

2.2.2 頻域特征提取

頻域特征需要從濾波后的多周期信號中提取，本文采用welch 法對小波分解濾波重構之后的脈搏信號序列進行功率譜估計，提取的特征包括諧波頻率、諧波幅值、諧波頻率差值df、譜能比、幅值差值、幅值之比和諧波面積之比，共33項特征[11]。

2.2.3 時-頻域特征提取

時-頻域聯(lián)合分析是信號處理分析中的重要部分，時-頻分析能將信號的局部特征表現(xiàn)出來，獲取信號在某一時刻的瞬時頻率特性，提取出不同于時域和頻域特征的重要特征信息。本文選用小波多尺度分析和經驗模態(tài)分解以及希爾伯特-黃變換三種常用的時-頻分析方法用于提取時-頻特征。

圖4 四種不同類型的脈象波形

信號f(t)在L2(R)空間上L2 范數(shù)的平方被定義為信號的能量，即若Dj(k)為信號f(t)經小波變換之后第j尺度上的第k個小波系數(shù)，則Ejk=|Dj(k)|2為信號在第j尺度k點上的小波能量，則Ej=為第j尺度N個采樣點的能量總和[13]。本文使用小波基Sym8 對多周期脈象序列進行10 層小波分解，獲取各尺度小波系數(shù)并依照上述定義計算了10層高頻細節(jié)分量以及第10 層低頻近似分量的小波系數(shù)能量值，將獲得的11個小波能量值進行歸一化即得到11個小波相關時-頻特征。

經驗模態(tài)分解（EMD）能夠根據(jù)信號本身的時間尺度來將其逐層分解成局部對稱的各個分量，它克服了傅里葉變換的局限性具有與小波變換相似的多分辨特性。經驗模態(tài)分解得到一系列固有模態(tài)函數(shù)（IMF）和一個殘余分量，能很好地保存信號本身的特性，避免信息丟失。希爾伯特-黃（Hilbert-Huang）變換是一種在非平穩(wěn)信號處理中應用廣泛的時-頻分析方法，克服了小波變換在時頻窗內依然具有平穩(wěn)性而導致在處理非平穩(wěn)信號時會產生諧波分量的缺點[14，15]。

本文對經過小波濾波后的多周期連續(xù)脈象序列進行經驗模態(tài)分解，篩選出5 個固有模態(tài)函數(shù)和一個殘余分量，如圖5 為一例經過小波分解濾波后的脈象信號的經驗模態(tài)分解結果。采用類似于上述小波能量特征的提取方法，提取5 個IMF 分量和1 個殘余分量的能量值，歸一化后得到6 個EMD 分解相關的時-頻能量特征。

對經過EMD 分解后的每一個IMF 分量進行希爾伯特-黃變換得到所有IMF 分量的Hilbert 譜，綜合所有IMF 分量的瞬時頻率并在時間軸上積分即可得到信號的Hilbert 邊際譜。Hilbert 邊際譜代表了各個頻率在全局幅值上的貢獻，表示統(tǒng)計意義上全部數(shù)據(jù)的累加幅度，如圖6 為上述同一例脈象信號的Hilbert 邊際譜。譜能比是一定頻率內能量Nj與總能量N的比值，由于脈象信號信息主要集中在0-20Hz，本文選取了 0.5、1、2、3、4、5Hz 內的能量與前 40Hz 內總能量的譜能比共6個特征。

圖5 脈象信號EMD分解結果

圖6 脈象信號Hilbert邊際譜

2.2.4 基于高斯混合模型的脈象信號特征提取

單周期脈搏波可以建模為高斯混合信號[5]，即為高斯信號，可以寫作gk(t)。ak、tk、σk分別是第k 個高斯信號函數(shù)的幅值、峰值位置和時間尺度因子，用三個高斯函數(shù)即可對單周期脈搏波形建模，可以得到共9個高斯混合模型參數(shù)特征。

綜上所述的特征提取方法，本文提取并用于評估分析的脈象特征如表2所示。

表2 提取脈象特征及其維數(shù)

3 基于隨機森林的脈象特征評估降維及結果

本文對脈象信號特征的研究和評估基于隨機森林算法（Random Forest,RF）。隨機森林是一種以決策樹為基學習器的集成學習算法[16]。首先對包含N 個樣本的訓練集進行M 輪自助采樣（有放回的重復獨立采樣），獲得M 個包含N 個訓練樣本的訓練集，基于這些訓練集訓練出M 棵未剪枝的決策樹，最終通過簡單投票法得出分類結果。

隨機森林在決策樹的訓練過程中引入了隨機屬性選擇。不同于傳統(tǒng)決策樹在選擇劃分屬性時每次選擇一個最優(yōu)屬性，在構建決策樹時會從該節(jié)點的屬性集合中隨機選取包含k 個屬性的子集，再從子集中選擇一個最優(yōu)屬性用于節(jié)點劃分，這種每次只需要考慮一個屬性子集的決策樹訓練方法使得RF 的訓練效率較高、容易實現(xiàn)且計算開銷小。

隨機森林在為節(jié)點選擇特征時根據(jù)Gini 增益最大化原理[17]，若父節(jié)點nf上的樣本被劃分到兩個子節(jié)點n1和n2中，則 Gini 增益最大化即是使 ΔIG=IG(nf)-p1*IG(n1)-p2*IG(n2)最大化，其中為節(jié)點 n 的 Gini 指數(shù)，pc為節(jié)點 n 中 c 類樣本所占比重，p1和p2則為從父節(jié)點nf劃分到兩個子節(jié)點n1和n2的樣本在父節(jié)點樣本中所占的比例。由于決策樹在生成節(jié)點時選擇能在該節(jié)點實現(xiàn)Gini 增益最大化的特征，因此特征的重要性可由節(jié)點樣本的劃分來體現(xiàn)，但RF 在訓練時引入數(shù)據(jù)樣本和輸入特征的雙重隨機性，可能導致區(qū)分度高的重要特征被用作劃分節(jié)點的次數(shù)比區(qū)分度低的特征少，所以不能簡單地用特征被用作劃分屬性的次數(shù)來衡量特征的重要性[18，19]。

特征si在第j 棵決策樹中作為節(jié)點劃分屬性的節(jié)點有N 個，則特征si在這棵決策樹上的特征重要性可表達為：

隨機森林中有J棵樹，特征si在整個RF中的重要性為：

3.1 基于隨機森林的特征評估可行性驗證

本文首先對基于隨機森林的特征重要性進行驗證，所選特征為使用前述特征提取方法所獲的時域、頻域、時-頻域和高斯混合信號特征，共93維。實驗中從樣本總量中分層抽樣其中的3/4 作為訓練集，剩余1/4 作為測試集，使得訓練集和測試集都保留了原有樣本總量中四種脈象數(shù)據(jù)類型的比例。為避免隨機森林算法中樣本抽取過程帶來的隨機因素影響，重復進行十次實驗。使用網格搜索法確定最佳的RF 參數(shù)，取每次測試集分類精度最高的一組作為特征選擇的依據(jù)，計算十組特征重要性結果的平均值并排序，實驗結果如表3。表中列出了按照特征重要性從高到低列出了部分特征的種類、名稱和特征重要性，其中前54維的累計重要性已達到95%。

根據(jù)表3 實驗結果，將排名前八十的特征按重要性降序分為四組，每組二十個特征。使用支持向量機（support vector machine, SVM）、反向傳播神經網絡（BP-NN）及隨機森林（random forest, RF）三種方法驗證特征重要性排序的正確性。本文的隨機森林算法采用Gini 指數(shù)作為判據(jù)，使用網格搜索法優(yōu)化算法的兩個關鍵參數(shù)：決策樹的數(shù)量和每棵決策樹的深度。SVM 模型性能主要由核函數(shù)、懲罰參數(shù)C 以及核函數(shù)參數(shù)γ決定，使用鏡像基函數(shù)作為核函數(shù)，并用網格搜索法優(yōu)化懲罰參數(shù)C 以及核函數(shù)參數(shù)γ。BP-NN 的網絡結構影響分類精度，目前并沒有明確的理論指導網絡層數(shù)與每層節(jié)點個數(shù)的選用。經過大量實驗，嘗試使用不同的網絡層數(shù)與節(jié)點個數(shù)，確定網絡結構為四層反向傳播神經網絡，具體為：兩個節(jié)點數(shù)為64 的全連接層，隨后一個Dropout層（系數(shù)為0.2）減少過擬合，最后使用softmax 層作為輸出層，其中優(yōu)化器為Adam。為避免訓練過程中的隨機性對實驗結論產生影響，每組特征組合均重復進行十次實驗，最后取十次實驗結果的平均值。

表3 基于隨機森林Gini指數(shù)的特征重要性

實驗結果如圖7 所示。從圖中可以看出，分類精度與隨機森林算法得到的特征重要性排序結果一致。三種算法的分類精度隨著特征重要性排名的下降而呈下降趨勢，特征的重要性排名越低，其分類精度越低。特征子集1 與特征子集4 的分類精度相差27.29%。

3.2 脈象特征降維及四類脈象分類實驗

在使用算法進行脈象分類前，通常會對較多的特征進行降維以提高算法的時間效率并通過去除可能導致錯分的冗余信息，減少過擬合提高分類準確率。常用的PCA 主成分分析將特征向量投影到新的坐標系中以降低維數(shù)，但經過變化后的特征向量無法用于研究原本特征在分類中的重要性，因此使用具體分類算法確定獲得最佳分類結果的維數(shù)有利于在獲取最佳分類精度的同時減少模型的復雜度，提高算法的運行效率。本研究采用序列前向選擇算法（Sequential Forward Selection,SFS），改變輸入特征的維度，計算相應RF、SVM以及BP-NN的分類精度。

結果如圖8 所示，初始時隨著重要性排名高的特征的加入，分類精度整體呈上升趨勢，隨后由于冗余特征的加入，分類精度隨后逐漸下降。結果證明本文的算法能夠有效減少脈象信號特征的數(shù)量，減少提取大量特征花費的時間，并能提高分類器的分類性能，三種算法達到最佳精度時，其特征數(shù)量都在15 維附近。

對于BP-NN 算法，其輸入特征分別在前15、16、18、19維獲得最佳精度，取最小維度15；對于RF算法，其最佳特征輸入維數(shù)為前13 維；對于SVM 算法；其最佳特征維數(shù)為前13 維。三種算法的特征類型如表4所示。

由表可以看出，在本文研究的平、實、弦、滑四種脈象分類中，時域特征對脈象信號的分類起主導作用，但同時排名靠前的時頻域特征和頻域特征也較為重要。

圖7 分類精度隨特征重要性排名的變化趨勢

圖8 分類精度隨特征數(shù)目變化趨勢

三種算法在最佳特征輸入維數(shù)與篩選前特征維數(shù)的分類精度對比如表5所示。表中精度提升的計算公式為。

表4 三種算法得到最佳精度時的特征組合

表5 三種算法特征篩選前后分類精度的提升

4 結語

人體脈象信號作為一種非平穩(wěn)信號，包含著與機體有關的多種生理信息，對于中醫(yī)脈診客觀化研究至關重要。本文對脈象信號常用的時域、頻域、時-頻域以及高斯特征進行了系統(tǒng)分析，并采用SVM、RF、BPNN 驗證基于Gini 指數(shù)的隨機森林算法可以得到特征重要性排名。采用序列前向選擇算法，改變輸入特征維度，當選用前15 個特征時，BP-NN 精度最高，相比篩選前分類精度提高12%，當選用前13 個特征時，SVM 和RF 精度均達到最高，相比降維前精度分別提高了13%和4.5%。本文進行分析的共175 例平、實、弦、滑四類脈象，時域特征在分類中占主導作用，同時基于頻域分析、小波分析和經驗模態(tài)分析的部分頻域和時-頻域特征也起到了重要作用。實驗結果表明本文算法可以根據(jù)特征重要性有效降低特征數(shù)量，確定不同分類算法適合的特征種類及維數(shù)并提高分類精度。