基于小波稀疏的心電信號壓縮感知?

劉繼忠嚴 旭鄧家誠華 晶

(1.南昌大學南昌市醫(yī)工結合技術研究重點實驗室，江西 南昌 330031；2.江西農業(yè)大學軟件學院，江西 南昌 33045)

心電信號(electrocardiograph，ECG)在心臟功能診斷上具有不可代替的作用[1]，但在經典奈奎斯特采樣理論下，采樣頻率必須高于待采樣信號最高頻率兩倍，才能無失真地恢復出原始信號。大量的采樣數據給可穿戴心電設備實時監(jiān)測與低功耗傳輸存儲帶來了不便。近年來出現的壓縮感知理論[2－3](compressed sensing，CS)打破了傳統采樣理論的限制，將采樣和壓縮融于一體，在信號稀疏域采集有效信息，為可穿戴心電設備的實時監(jiān)測和低功耗傳輸提供了新的解決思路。

在CS理論框架下，實時監(jiān)護系統可在心電信號稀疏域采樣、達到低數據量的目的，后續(xù)可通過重構算法恢復原信號或者直接在壓縮域下完成識別以進一步達到低功耗目標[4－5]。文獻[6]中作者對比分析多種離散小波基的重構性能，對心電信號壓縮比可達70%，不足在于忽略心電信號投影到不同小波基的稀疏度問題，重構均方根差百分比誤差(percentage root-mean-squared difference，PRD)高于10%。文獻[7]中基于時空稀疏模型對心電信號壓縮采樣，利用字典學習構建符合心電信號特點的稀疏矩陣，在壓縮比為70%時、實現重構誤差9.03%，實現高壓縮比和較低重構誤差，缺點在于未考慮觀測矩陣對重構精度的影響。文獻[8]基于心電信號波形特點，提出塊稀疏貝葉斯學習算法，但使用的重構算法單一，信號恢復效果仍不夠理想。

因此，文中依據CS理論，分別從投影稀疏度角度探討小波域下心電信號的最優(yōu)稀疏表示，以及觀測矩陣和算法對重構誤差的影響，并以MIT-BIH數據庫中心電信號作為原始信號，綜合研究在高壓縮比下實現低誤差重構，達到可穿戴實時健康監(jiān)護系統低數據量采集傳輸和終端高精度重構的目的。

1 心電信號稀疏表示

1.1 信號稀疏表示

CS要求信號本身稀疏或在某個變換域稀疏，而自然界大多數信號本身并非稀疏，心電信號在時域也非稀疏，因此需要設計合適的稀疏變換矩陣，得到心電信號投影到稀疏矩陣上的稀疏表示。考慮到可穿戴心電監(jiān)護設備嵌入式硬件系統的實現和小波在心電信號壓縮[9]、去噪中的良好性能[10]，將重點探討心電信號離散小波變換(discrete wavelet transform，DWT)下的稀疏表示[11]，取daubechies(dbN)、coiflet，(coifN)、symlet(symN)、biorthogonal(biorNr.Nd)和reverseBior(rbioNr.Nd)等可離散變換的小波構建變換矩陣，并取MIT-BIH數據庫第100號心電信號前1 024個點作為實驗數據，以同源小波濾波去噪后的信號作為實驗信號X1024(X∈RN×1，N＝1 024)。同源小波濾波即濾波小波與構建稀疏矩陣小波采用同一母小波函數[12]。

對心電信號X作離散序列Mallat分解，構建離散小波變換矩陣即CS中稀疏矩陣如式(1)、式(2):

式中:G(n)、H(n)為所選取小波函數對應的低通和高通濾波器系數序列，G(n)，H(n)∈RL×1。j為分解層數，j＝0時，X0(n)為實驗心電信號X1024，卷積后G1(n)∈R(N＋L－1)×N，將前L/2－1行和后L/2－1行去除，從剩余N＋1行中挑選偶數行共N/2行，H1(n)同理挑選N/2行，疊加組成一層分解的小波變換矩陣Ψ∈RN×N。從H1(n)中挑選出的N/2行組成低頻部分X1(n)，依式(1)、式(2)繼續(xù)迭代可得多層分解小波變換矩陣。則X稀疏表示如式(3)。

式中:X為心電信號，定義Ψ為稀疏矩陣，S為心電信號稀疏表示，K為稀疏度，S∈RN×1，S中只有K個分量大于閾值γ，γ→0＋。心電信號投影到不同的稀疏矩陣上，信號稀疏度越低，需采樣的數據越少，有利于數據存儲和傳輸。

1.2 稀疏表示分析

不同分解層數和不同小波母函數將對應不同的稀疏矩陣Ψ，依據式(3)得到的稀疏表達S會有差異。取閾值γ＝0.05，將心電信號投影到Daubechies小波、Coiflet小波、Symlet小波、Biorthogonal小波和ReverseBior小波5種小波不同分解層數對應的稀疏矩陣上，得到稀疏度匯總結果如圖1所示。

由圖1可知，隨著分解層數的增加，稀疏度呈降低趨勢，在分解層數達到第7層時各離散小波心電信號的稀疏變換分解基本達到了穩(wěn)定值，而分解到第9層時全部達到了穩(wěn)定值。圖1(a)中，Db4小波在第6層時其稀疏變換矩陣對應的稀疏度就達到了穩(wěn)定，且稀疏度在Daubechies小波中最低，對應稀疏度K＝98；圖1(b)中，Coif4小波在第5層時就達到了穩(wěn)定，且稀疏度在Coiflet小波中最低，對應K＝104；圖1(c)中，Sym5小波在第6層時就達到了穩(wěn)定，且稀疏度在Symlet小波中最低，對應K＝94；圖1(d)中，Bior5.5小波在第6層時達到了穩(wěn)定，且稀疏度在Biorthogonal小波中最低，對應K＝100；圖1(e)中，Rbio1.3小波在第7層時達到了穩(wěn)定，且稀疏度在ReverseBio小波中最低，對應K＝89。同時，在相同條件下，X1024經過離散正弦基(discrete sine transform，DST)稀疏表示后，K＝252，經過離散余弦基(discrete cosine transform，DCT)稀疏表示后，K＝245。

圖1 心電信號不同離散小波域稀疏度

因此，從稀疏度角度而言，可從對應離散小波系中取Db4小波、Coif4小波、Sym5小波、Bior5.5小波、Rbio1.3小波用于構建稀疏變換矩陣，依次可得到對應小波域中心電信號最低稀疏度的稀疏表示，而心電信號投影到DST、DCT后K值過高，不宜作為稀疏矩陣。

2 觀測矩陣的構建

2.1 觀測矩陣概念

觀測矩陣Φ(Φ∈RM×N，M?N)起到將N維信號降維到M維的作用，M值的確定如式(4)。通過觀測矩陣Φ將原信號X1024壓縮采樣得到測量值Y，如式(5)。由式(3)、式(5)，可得式(6)。而在后續(xù)重構階段，稀疏矩陣Ψ，觀測矩陣Φ和測量值Y作為先驗條件，通過重構算法從M維測量值Y恢復出N維稀疏近似，進而恢復整段心電信號X1024。

觀測矩陣Φ同時完成數據采樣和壓縮，其須滿足約束等距特性(restricted isometry property，RIP)，在重建階段能夠高概率近乎完整恢復出原始信號[13]。RIP準則如式(7)所示:

式中:ζS為有限等距常數[14](restricted isometry constant，RIC)，要求0≤ζS≤1。

2.2 觀測矩陣分析

由于RIP準則中有限等距常數較難計算驗證，Tao等人證明了伯努利矩陣(Bernoulli matrix)、循環(huán)矩陣(Circulant matrix)、高斯隨機矩陣(Gaussian matrix)、部分哈達馬矩陣(PartHadamard matrix)和托普利茲矩陣(Toeplitz matrix)等大概率滿足RIP條件[15－17]。

為節(jié)約計算分析時間，文中進一步取心電數據X1024中的前256位數據作為實驗信號，記為X256。依據2中稀疏表達結論，以Db4小波、Coif4小波、Sym5小波、Bior5.5小波、Rbio1.3小波構建稀疏矩陣，選取前述高斯隨機矩陣、伯努利矩陣、循環(huán)矩陣、部分哈達馬矩陣、托普利茲矩陣作為觀測矩陣研究對象，并采用經典正交匹配追蹤算法(orthogonal matching pursuit，OMP)完成信號重構(4中將進一步討論信號重構)。假設重構信號256與實驗信號X256之間殘差小于0.1為重構成功，定義重構成功率如式(8)所示:

確定稀疏矩陣和重構方法后，在不同的觀測矩陣下逐次增大M(與壓縮比有關)，直到M＝＝N。文中取M為一定值時的重構實驗次數為100，即重復重構一百次并記錄成功次數得到重構成功率。

分別以Db4小波、Coif4小波、Sym5小波、Bior5.5小波、Rbio1.3小波構建稀疏矩陣，繪制得到不同觀測矩陣下觀測值M與重構成功率的關系曲線，如圖2所示。

圖2 不同稀疏矩陣和觀測矩陣下的心電信號重構成功概率

由圖2可知，重構成功率隨M的增高(即壓縮比的降低)而顯著提高，對于不同的小波稀疏，高斯觀測矩陣和伯努利觀測矩陣的重構性能明顯優(yōu)于其他矩陣，而伯努利觀測矩陣比高斯矩陣稍提前實現重構成功率達到100%，且矩陣元素為更易于硬件實現的±1。當測量矩陣為伯努利矩陣時，相比較而言，Db4小波、Rbio1.3小波在較低觀測值(M約80，Db4小波M＝81，Rbio1.3小波M＝80)即較高的壓縮比時，可達到100%成功率；而Sym5小波次之，M＝100時達到100%成功率。即以Db4小波、Rbio1.3小波構建稀疏矩陣能以更少的測量值(更高的壓縮比)實現高重構成功率。

因此，綜合上述分析并考慮到嵌入式硬件系統的實現問題，稀疏矩陣以Bb4小波和Rbio1.3小波構建，觀測矩陣選用伯努利矩陣，以實現少測量值高概率重構心電信號的目的。

3 心電信號的重構

3.1 CS信號重構方法

要實現心臟功能的觀察診斷，需要在終端實現心電信號的重構。重構即是以測量值向量Y，借助稀疏矩陣、觀測矩陣和重構算法恢復N維心電信號的過程。

重構過程為式(6)的逆向求解過程，由于M?N，求解方程為欠定方程，但S是K－稀疏的，當觀測矩陣Φ滿足RIP準則時，求解可轉換為?0范數約束優(yōu)化問題，如式(9)。對應2中尋找最稀疏的分解表達，稀疏度越低，重建效果越好，逆向求解得到稀疏表達的近似值＾S，然后依據式(3)得到重構信號。

壓縮感知重構方法有貪婪算法和凸優(yōu)化算法，貪婪算法主要包含OMP算法及其改進型；凸優(yōu)化算法以基追蹤法(basis pursuit，BP)為代表。本部分將探討在前述稀疏變換和觀測矩陣優(yōu)化的基礎上適合心電信號的恢復重構方法。

3.2 基本型OMP算法

步驟如下:

①初始化殘差。

②從感知矩陣θ(θ＝{θ1，θ2，…，θN})中選出與R i最相關的原子θi滿足式，r0為感知矩陣列索引。信號Y被分解成矩陣投影分量和殘差R1。

③將選出的最相關原子θr0添加到矩陣Λi中。

④將R i投影到矩陣Λi所張成空間上得到投影值Pi。

⑤更新殘差。

循環(huán)步驟②～⑤直到i＝k時停止。用最小二乘法得到K－稀疏近似值，如式

3.3 凸優(yōu)化算法

基本原理是將求解?0范數問題放寬到?1范數凸優(yōu)化問題，如式(17)所示，使整個重構能用線性規(guī)劃來解決，BP優(yōu)化可借助L1＿magic工具箱實現。

3.4 改進型OMP算法及凸優(yōu)化類算法重構評估

結合正交匹配追蹤算法的特點，對正交匹配追蹤算法提出了不同改進[18]:

①廣義正交匹配追蹤(generalized OMP，GOMP):在式(10)中，每次迭代不再選出單個原子，而是按內積大小選出前t個，更新Λi，得GOMP算法；

②稀疏度自適應匹配追蹤算法(sparsity adaptive MP，SAMP):在式(11)中，算法不再運行K次，而是設置殘差閾值ε，如果新殘差Ri小于ε，則運行結束，得SAMP算法；

③分段正交匹配追蹤算法(stagewise OMP，StOMP):在式(10)中，選擇原子時設置門限閾值，內積大于閾值，將原子選出并更新Λi。同時設置循環(huán)次數，算法不機械式循環(huán)K次，得StOMP算法；

④正則化正交匹配追蹤算法(regularized OMP，ROMP):在式(10)中，按內積選出前K個原子，從K個原子中選出與殘差內積絕對值的最大值不能比最小值大兩倍以上且能量最大的一組原子，循環(huán)K次，得ROMP算法；

⑤壓縮采樣匹配追蹤算法(compressive sampling MP，CoSaMP):在ROMP基礎上，不斷更新選出的原子，上一次迭代選出原子下次可能會被頂替，得到CoSaMP算法。

為評價心電信號壓縮采樣與重構的整體性能，使用壓縮比(compression ratio，CR)表征壓縮程度，CR越高壓縮性越好，如式(18)所示。采用均方根差百分比表征重構信號與原始信號之間誤差，PRD越小表明重構誤差越小，如式(19)所示。

依據前述分析，實驗中稀疏矩陣以db4小波和rbio1.3小波構建，觀測矩陣采用伯努利矩陣。由式(4)可知，M與C和K取值相關，由2中結果取db4小波和rbio1.3穩(wěn)定分解后稀疏值K分別為98和89，且經實驗分析此情況下C≥2.5時，M的取值使PRD穩(wěn)定。因此，我們取C＝2.5對應的M值進行運算分析，在Win10 Intel Pentium 3556U 1.7 GHz和MATLAB 2016平臺上分別運行OMP及其改進型算法、及求解?1范數的BP凸優(yōu)化方法，得到重構心電信號＾X1024，各重構算法隨機運行10次，得到10次運行平均參數結果如表1所示。

由表1可知，在C固定(二者CR指數接近)的情況下，Db4小波稀疏的重構效果在7種算法中都優(yōu)于Rbio1.3小波，而且在GOMP、SAMP的重構實驗中，從PRD指數即重構誤差上看，以Db4小波稀疏遠優(yōu)于Rbio1.3小波。

表1 不同算法下心電信號壓縮采樣與恢復重構性能

從總的重構精度角度分析，ROMP算法在心電信號兩種小波稀疏的CS重構中誤差較大，同時，對于Rbio1.3小波稀疏，其余算法的重構誤差相近(在同一個數量級上)，對于Db4小波稀疏，GOMP、SAMP算法的重構誤差較低，OMP、StOMP、CoSaMP和BP稍高。從重構耗時角度，BP、SAMP算法重構耗時較長，CoSaMP耗時次之，OMP、GOMP、StOMP耗時較短，而ROMP耗時最少。

綜合考慮壓縮采樣效果、重構精度及算法重構計算效率，基于Db4小波稀疏、伯努利觀測矩陣、GOMP重構算法的心電信號壓縮感知采樣方法為最優(yōu)，可以實現高的壓縮效率、重構精度和重構效率。

3.5 壓縮采樣測量值及壓縮比

為進一步分析文中提出的基于Db4小波稀疏、伯努利觀測矩陣采樣、GOMP算法重構的心電信號壓縮感知方法，從而滿足實際應用中心電信號數據壓縮率和重構精度的需求，考慮壓縮感知采樣的CR指數與PRD指數，給出上述采樣方法下M從300～480變化的X1024壓縮采樣及重構的M－PRD關系分析，結果如圖3所示。

圖3 M－PRD關系圖

由圖3可知，隨著測量值M取值的增加，重構誤差PRD呈下降趨勢，在M＝315時，重構誤差PRD約為0.996 2%，M＝450時，PRD僅為0.001 4%。以取M＝320為例，CR＝68.75%(保留的有效數據約為30%)，重構誤差PRD＝0.914 2%，原信號與重構信號的對比如圖4，由圖4可知重構信號基本完整覆蓋原信號。

圖4 重構對比圖

為進一步驗證文中提出方法的有效性，隨機選取MIT-BIH數據庫第105號、119號、122號、201號、220號五組心電信號前1 024點作為新實驗數據，取M＝320，CR＝68.7%，采用Db4小波濾波及稀疏、伯努利觀測矩陣、GOMP重構算法，重復上述稀疏、觀測、重構步驟10次并記錄PRD指數以及10次均值及其標準差，得到實驗結果如表2所示。

表2 不同心電信號壓縮傳感采樣與重構實驗結果

由表2可知，對于1 024點的不同組心電信號，在測量數為320(壓縮比為68.75%)的情況下，各信號的PRD誤差皆小于或約等于1%，重構耗時在0.2 s左右，驗證了所選擇方案的可行性與穩(wěn)定性，可以較好的實現可穿戴心電監(jiān)控設備小數據量低功耗傳輸(保留約30%數據)和高精度重構(重構PRD誤差小于或約等于1%)，并可根據實際要求調整CR和PRD(或者硬件性能許可情況下設計自適應信號稀疏特性的CR)，具有較強的實際應用價值。

4 結束語

為可穿戴心電信號實時監(jiān)測設備實現小數據量、低功耗傳輸和高精度重構，文中以MIT-BIH心電信號數據庫樣本為例，探討了離散小波稀疏域與不同觀測矩陣、不同重構算法下的心電信號稀疏變換、觀測矩陣壓縮采樣、高精度重構等壓縮感知問題，進行了綜合分析研究和實驗驗證，提出了基于db4小波濾波稀疏、伯努利觀測矩陣采樣、GOMP重構算法的心電信號壓縮感知采樣方法。結果表明:

①同源濾波下，離散小波系中Db4小波、Coif4小波、Sym5小波、Bior5.5小波、Rbio1.3是同類小波中獲得最低稀疏度的小波稀疏基。

②對于不同小波稀疏，高斯觀測矩陣和伯努利觀測矩陣的重構性能明顯優(yōu)于其他矩陣，而在實驗條件下伯努利觀測矩陣比高斯矩陣稍提前實現重構成功率達到100%，以Db4小波、Rbio1.3小波構建稀疏矩陣能以更小的測量值(即更高的壓縮比)實現高重構成功率，且更易于硬件實現。

③在文中7種重構算法中，Db4小波稀疏的重構效果都優(yōu)于Rbio1.3小波；而以Db4小波稀疏的重構算法中GOMP、SAMP重構算法精度較高，OMP、GOMP、StOMP耗時較少，ROMP耗時最少，綜合提出了基于Db4小波稀疏、伯努利觀測矩陣、GOMP重構算法的心電信號壓縮感知采樣與重構方法。

④CR＝68.75%(約保留30%有效稀疏數據)情況下，心電信號重構PRD誤差小于或約等于1%，隨機抽取數據庫其他心電信號實驗驗證了方法的可行性和穩(wěn)定性，為可穿戴心電設備實時低功耗監(jiān)測傳輸實現提供了支撐。