車波，賁鴻偉，朱霖霖，劉磊，鄧林紅 △

(1.常州大學(xué)生物醫(yī)學(xué)工程與健康科學(xué)研究院,常州 213164;2.常州大學(xué)信息科學(xué)與工程學(xué)院,常州 213164)

1 引 言

小氣道功能檢查在氣道阻塞性疾病中具有重要的評(píng)估意義[1-3]。在現(xiàn)有的肺功能檢查方法中，強(qiáng)迫振蕩技術(shù)(forced oscillation technique, FOT)是一種無(wú)創(chuàng)、快速和配合要求少的檢查方法[4-5]。FOT自1956年由DouBois等首次引入肺功能檢查以來(lái)，目前已經(jīng)廣泛應(yīng)用于肺功能輔助診斷[6-7]。研究表明，低頻的振蕩激勵(lì)信號(hào)能更多地反映小氣道組織的流變學(xué)特性及阻塞情況[8-10]。不過(guò)，由于該低頻段接近人體呼吸頻段 (0.1～1.1 Hz)，易發(fā)生頻帶的混疊干擾，并且該非線性干擾很難用現(xiàn)有FOT技術(shù)來(lái)分離和消除[11-12]。因此，有必要在低頻激勵(lì)的FOT檢測(cè)中，引入一種能有效處理非線性、非平穩(wěn)信號(hào)的分解和降噪方法。

目前非平穩(wěn)信號(hào)的自適應(yīng)分解方法， 主要包括經(jīng)驗(yàn)?zāi)B(tài)分解(empirical mode decomposition，EMD)及改進(jìn)方法[13-14]、變分模態(tài)分解(variational mode decomposition，VMD)[15-16]等。其中，EMD及改進(jìn)方法的數(shù)學(xué)理論基礎(chǔ)薄弱且計(jì)算量較大。相比之下，VMD既有嚴(yán)格的數(shù)學(xué)理論支撐，其構(gòu)建的變分優(yōu)化模型相對(duì)簡(jiǎn)潔。最近，Maji等[17]對(duì)比EMD和VMD在心電噪聲中的分離效果，表明VMD能更有效地去除高頻噪聲。Nazari等[18]則利用VMD成功地從心電圖中提取到呼吸率信號(hào)。但目前對(duì)低頻激勵(lì)條件下呼吸阻抗信號(hào)非線性分解和降噪的研究還較少，僅見Copot等[19]采用EMD方法減小呼吸阻抗偏移的報(bào)道，也并未比較不同的分析方法。

鑒于此，本研究引入VMD方法用于低頻激勵(lì)下的呼吸壓力和流量信號(hào)的降噪，并對(duì)比分析VMD與EMD在阻抗信號(hào)分解和降噪中的優(yōu)勢(shì)，進(jìn)而探討VMD降噪的參數(shù)優(yōu)化和低頻振蕩激勵(lì)等問(wèn)題。

2 實(shí)驗(yàn)方法

2.1 呼吸信號(hào)采集與處理系統(tǒng)

本研究采用實(shí)驗(yàn)室自主設(shè)計(jì)制作的FOT系統(tǒng)來(lái)測(cè)量呼吸壓力和流量信號(hào)，其工作原理見圖1。

圖1 FOT系統(tǒng)示意圖

該系統(tǒng)主要由激勵(lì)源、傳感器和上位機(jī)處理三大模塊構(gòu)成。激勵(lì)源模塊采用直線音圈電機(jī)激勵(lì)源，可產(chǎn)生1～5 Hz頻率、0.2 ～ 0.98 kPa壓力的單正弦和方波脈沖激勵(lì)U(t)。在傳感器和上位機(jī)處理模塊中，利用LabVIEW設(shè)計(jì)的界面可實(shí)時(shí)處理來(lái)自傳感器采集的壓力P(t)、流量Q(t)信號(hào)。系統(tǒng)各參數(shù)按照之前發(fā)表的數(shù)值進(jìn)行設(shè)置[20]。

圖2 FOT系統(tǒng)辨識(shí)原理圖

FOT辨識(shí)系統(tǒng)見圖2，呼吸總阻抗Zrs的計(jì)算公式：

(1)

P(ω)、Q(ω)分別為頻域上的壓力和流量，實(shí)部Rrs為呼吸阻力，Xrs為呼吸電抗。

2.2 實(shí)驗(yàn)數(shù)據(jù)采集對(duì)象

本研究數(shù)據(jù)均來(lái)自實(shí)驗(yàn)室志愿受試者，共測(cè)試了3名健康成年男性，具體參數(shù)見表1。實(shí)驗(yàn)參照歐洲呼吸學(xué)會(huì)(ERS)的推薦標(biāo)準(zhǔn)，記錄受試者自主平靜呼吸下三個(gè)呼吸周期(約60 s)的信號(hào)，且重復(fù)5次。取其中波形穩(wěn)定、無(wú)較大呼吸偽影的記錄作為最終實(shí)驗(yàn)數(shù)據(jù)。

表1 采集對(duì)象實(shí)驗(yàn)數(shù)據(jù)

2.3 VMD信號(hào)分解方法

本研究中的VMD方法是基于復(fù)雜信號(hào)頻域稀疏性的假設(shè)，通過(guò)構(gòu)造約束性變分模型，迭代搜索最優(yōu)解，實(shí)現(xiàn)復(fù)雜信號(hào)的自適應(yīng)頻域窄帶分量的分離。

2.3.1VMD算法的分解流程 VMD算法對(duì)混合信號(hào)的分解流程見圖3。變量和參數(shù)為呼吸壓力信號(hào)x(t)及各模態(tài)u(t)，分量中心頻率ω，懲罰因子α和Lagrange乘子λ，模態(tài)數(shù)K，收斂閾值ε。

其中，最優(yōu)解搜索的關(guān)鍵步驟：

圖3 混合信號(hào)的VMD分解流程圖

(2)

(3)

步驟2：對(duì)所有ω≥0，進(jìn)行雙重提升。

(4)

本研究利用MATLAB 編寫myVMD函數(shù)，[u, u_fres, omega]= myVMD(signal, alpha, tau,K, DC, init, maxIters, tol)，參數(shù)說(shuō)明見表2。

表2 myVMD的參數(shù)說(shuō)明

2.3.2模態(tài)數(shù)K的確定方法 為有效抑制VMD分解的模態(tài)混疊，需選取恰當(dāng)?shù)哪B(tài)數(shù)K。在合適的K值范圍內(nèi)，其瞬時(shí)頻率的均值變化較平緩具有物理意義；當(dāng)K較大時(shí)，分量瞬時(shí)頻率的均值曲線會(huì)出現(xiàn)顯著下跌，故可通過(guò)求解分量瞬時(shí)頻率均值曲線曲率的方法來(lái)確定合適的K值。

此外，本研究采用的計(jì)算實(shí)驗(yàn)環(huán)境為：Windows7 64bit，MATLAB R2017b。myVMD方法的輸入?yún)?shù)設(shè)置為：alpha=1 200, tau=0,K=3, DC=1, init=1, maxIters=500, tol=1e-5。其中，alpha取采樣頻率fs(600 Hz)的兩倍，tau在當(dāng)前強(qiáng)噪聲背景下設(shè)為0。

2.4 EMD方法說(shuō)明

EMD方法是通過(guò)迭代篩分的方式，將復(fù)雜信號(hào)x(t)自適應(yīng)分解成多個(gè)單分量(本征模態(tài)函數(shù)，imfi(t))和殘余分量rn(t)：

r1(t)=x(t)-imf1(t)

(5)

(6)

在方法對(duì)比中，由于基于經(jīng)驗(yàn)框架的EMD及改進(jìn)方法與基于變分框架的VMD不同。因此，為對(duì)比這兩類方法對(duì)呼吸信號(hào)的分解效果，選取了EMD和VMD方法的處理效果作比較分析。

3 實(shí)驗(yàn)結(jié)果

圖4列出了一段呼吸壓力信號(hào)經(jīng)VMD分解后的各分量瞬時(shí)頻率的均值曲線，當(dāng)K增加到4時(shí)，曲線明顯下降，故后續(xù)VMD分解計(jì)算中取K=3。

圖4 分量瞬時(shí)頻率的均值變化曲線

圖5為受試者X1的一段呼吸壓力信號(hào)及其經(jīng)EMD分解結(jié)果，按頻率降序排列依次為分量IMF1-8和殘余分量Res。

圖6為同一呼吸壓力信號(hào)經(jīng)過(guò)VMD分解的結(jié)果，其中分量BLIMF3已經(jīng)能很好地反映混合信號(hào)的較小波動(dòng)趨勢(shì)。

圖5 EMD分解壓力信號(hào)的IMF 1-8和Res分量

圖6 VMD分解壓力信號(hào)的BLIMF 1-3分量

對(duì)比圖5和圖6的結(jié)果可知， VMD方法在K=3時(shí)，可以很好地分解低頻部分，需分解的分量數(shù)也較少。EMD方法由于分量個(gè)數(shù)不可控，因而計(jì)算量較大。盡管將EMD多個(gè)分量疊加也能實(shí)現(xiàn)降噪效果，但需要疊加的分量不能預(yù)知。因此，VMD方法對(duì)本研究采集的呼吸信號(hào)中低頻信號(hào)的分解具有明顯優(yōu)勢(shì)。

VMD方法對(duì)呼吸壓力信號(hào)降噪的結(jié)果見圖7，可見VMD方法在強(qiáng)噪聲背景(灰色)下，通過(guò)分離出中高頻信號(hào)，得到清晰的低頻信號(hào)(紅色為降噪后流量，藍(lán)色為降噪后壓力)，較好地保留了低頻信號(hào)的局部非線性特征。

圖7 壓力(藍(lán)線)、流量(紅線)經(jīng)VMD降噪后波形

圖8為呼吸總阻抗功率譜密度的時(shí)頻分布，顯示了呼吸總阻抗在低頻段的非平穩(wěn)性。其功率譜密度在0～15 Hz范圍內(nèi)隨時(shí)間變化的波動(dòng)較大，尤其在5 Hz及以下頻段表現(xiàn)出較高的非平穩(wěn)性。這些為分析低頻段阻抗的頻率依賴性提供了依據(jù)。

圖8 呼吸總阻抗Zrs的功率譜密度

圖9為經(jīng)VMD處理計(jì)算5個(gè)激勵(lì)頻率下的25個(gè)阻抗值分布，由實(shí)軸Rrs和虛軸Xrs構(gòu)成分布的復(fù)平面，兩軸上的數(shù)值近似呈高斯分布。線性回歸(紅線)分析了阻抗值集中在Rrs[0.25, 0.4]及Xrs[-0.3, -0.1]范圍內(nèi)。圖9中阻抗值過(guò)小的點(diǎn)(綠圈內(nèi))，可能是由于測(cè)量中的泄漏造成，在計(jì)算中可舍棄，以減小該人工呼吸偽影的影響。

圖9 Zrs在Rrs和Xrs復(fù)平面的散點(diǎn)分布

圖10為VMD和EMD方法對(duì)呼吸阻抗估計(jì)的均值曲線。Rrs(上欄)和Xrs(下欄)結(jié)果顯示，兩種方法的處理結(jié)果在1 Hz和5 Hz端點(diǎn)處的方差相對(duì)較大，但均小于 0.1 kPa·s/L 。Rrs和Xrs阻抗在1～5 Hz具有一定的頻率依賴性，且變化趨勢(shì)連續(xù)。相比EMD處理的結(jié)果，經(jīng)VMD處理后的阻抗曲線變化更為平緩。阻抗頻譜結(jié)果也表明，VMD方法有助于在低頻激勵(lì)條件下準(zhǔn)確估計(jì)Rrs和Xrs。由Xrs曲線可初步預(yù)測(cè)其共振頻率Fres(過(guò)零點(diǎn)的頻率)在6 Hz附近，符合健康成年人的Fres低于10 Hz的參考標(biāo)準(zhǔn)。

圖10 EMD(藍(lán)線)和VMD(紅線)方法的阻抗頻譜比較

4 結(jié)果與討論

本研究的VMD方法在低頻振蕩激勵(lì)條件下對(duì)呼吸阻抗測(cè)量的分解和降噪應(yīng)用，實(shí)現(xiàn)了對(duì)呼吸壓力和流量信號(hào)的自適應(yīng)分解。阻抗分析結(jié)果表明：VMD方法對(duì)低頻段阻抗的分解效果比EMD方法更高效、穩(wěn)定，這為后面準(zhǔn)確估計(jì)低頻段下的小氣道阻抗提供了分析基礎(chǔ)。

此外，有兩個(gè)問(wèn)題值得討論和研究。第一，VMD中參數(shù)K和α對(duì)分解效果的影響較大，其值應(yīng)由信號(hào)本身來(lái)估計(jì)，以實(shí)現(xiàn)數(shù)據(jù)自驅(qū)動(dòng)分解。本研究分別利用了分量瞬時(shí)頻率的均值和采樣率來(lái)估計(jì)K、α，但對(duì)于兩者的組合最優(yōu)化問(wèn)題還需進(jìn)一步研究。第二，考慮氣道阻抗的頻率依賴性，即在0.1～1 Hz頻段激勵(lì)下的阻抗檢測(cè)問(wèn)題。由于該低頻段更接近正常人體呼吸頻率，可更準(zhǔn)確反映小氣道的粘彈性變化。若能在更低頻激勵(lì)條件下，結(jié)合阻抗的力學(xué)參數(shù)模型和VMD等自適應(yīng)分解方法，對(duì)于準(zhǔn)確評(píng)估低頻激勵(lì)下的呼吸阻抗具有重要意義。