基于貝葉斯濾波方法的心臟運(yùn)動(dòng)信號(hào)實(shí)時(shí)濾波算法

梁 帆，楊莉莉，崔世鋼，吳興利

（天津職業(yè)技術(shù)師范大學(xué)天津市信息傳感與智能控制重點(diǎn)實(shí)驗(yàn)室，天津 300222）

基于貝葉斯濾波方法的心臟運(yùn)動(dòng)信號(hào)實(shí)時(shí)濾波算法

梁 帆，楊莉莉，崔世鋼，吳興利

（天津職業(yè)技術(shù)師范大學(xué)天津市信息傳感與智能控制重點(diǎn)實(shí)驗(yàn)室，天津 300222）

提出一種基于貝葉斯濾波遞推狀態(tài)估計(jì)的心臟運(yùn)動(dòng)測(cè)量信號(hào)實(shí)時(shí)濾波算法，提高心臟數(shù)據(jù)的測(cè)量精度，該方法建立了心臟運(yùn)動(dòng)模型與傳感器誤差模型，采用了擴(kuò)展卡爾曼濾波與粒子濾波兩種貝葉斯濾波方法的實(shí)現(xiàn)形式。實(shí)驗(yàn)結(jié)果表明，提出的算法有效地還原了心臟運(yùn)動(dòng)的真實(shí)值，并且滿足了跟蹤系統(tǒng)實(shí)時(shí)數(shù)據(jù)處理的要求，為整個(gè)系統(tǒng)精度的提高提供了保證。

手術(shù)輔助機(jī)器人；實(shí)時(shí)濾波；擴(kuò)展卡爾曼濾波器；粒子濾波器

在機(jī)器人輔助非體外循環(huán)旁路移植（coronary artery bypass grafting，CABG）手術(shù)中，手術(shù)輔助機(jī)器人的主要作用是實(shí)時(shí)跟蹤心臟運(yùn)動(dòng)，保持心臟表面手術(shù)點(diǎn)與機(jī)器人操作端點(diǎn)之間的相對(duì)靜止，幫助醫(yī)生完成對(duì)于人手難以實(shí)現(xiàn)的快速、高幅度心臟運(yùn)動(dòng)的跟蹤任務(wù)[1]。為了實(shí)現(xiàn)精確預(yù)測(cè)與跟蹤控制，低噪和高質(zhì)的心臟運(yùn)動(dòng)測(cè)量數(shù)據(jù)是手術(shù)輔助機(jī)器人控制系統(tǒng)的前提基礎(chǔ)和重要保證。在以往的心臟手術(shù)輔助機(jī)器人相關(guān)研究中，大多數(shù)情況下研究者首先對(duì)采集到的心臟運(yùn)動(dòng)數(shù)據(jù)做離線處理，然后應(yīng)用處理后的數(shù)據(jù)作為跟蹤控制算法的參考數(shù)據(jù)，用這些數(shù)據(jù)來進(jìn)一步驗(yàn)證運(yùn)動(dòng)跟蹤算法[2-3]。但是在真實(shí)的心臟手術(shù)過程中，測(cè)量數(shù)據(jù)必須以實(shí)時(shí)在線的方式進(jìn)行處理，并傳送給后端的控制器作為參考輸入。本文提出一種在噪聲干擾環(huán)境下對(duì)心臟運(yùn)動(dòng)信號(hào)實(shí)時(shí)濾波的可靠方法。

心臟運(yùn)動(dòng)信號(hào)的采集通過超聲波距離測(cè)量系統(tǒng)。在真實(shí)環(huán)境中，超聲傳感系統(tǒng)主要的誤差產(chǎn)生源是水晶振子空間尺寸準(zhǔn)確度以及超聲波回聲干擾。如果直接采用超聲傳感系統(tǒng)測(cè)量值，計(jì)算出的信號(hào)往往由于噪聲污染的不確定性而產(chǎn)生嚴(yán)重誤差。因?yàn)樵趯?shí)際的心臟手術(shù)中，心臟運(yùn)動(dòng)過程中超聲波傳遞可能會(huì)受到中間物的阻攔而使測(cè)量位置產(chǎn)生偏差，從而直接影響到后續(xù)跟蹤控制的參考信號(hào)質(zhì)量。本文中用基于貝葉斯濾波器[4]的狀態(tài)估計(jì)方法對(duì)傳感器測(cè)量數(shù)據(jù)進(jìn)行信號(hào)處理，利用測(cè)量值的統(tǒng)計(jì)特征從觀測(cè)數(shù)據(jù)中恢復(fù)測(cè)量值的真實(shí)狀態(tài)。其中將心臟運(yùn)動(dòng)信號(hào)的真實(shí)值看作狀態(tài)變量，通過與心臟運(yùn)動(dòng)模型相關(guān)的相鄰時(shí)刻產(chǎn)生狀態(tài)轉(zhuǎn)移概率和包含測(cè)量誤差統(tǒng)計(jì)特征的測(cè)量概率的每個(gè)時(shí)刻的更新，遞推地估計(jì)出測(cè)量真值，力求有效降噪并減低測(cè)量過程中不確定度帶來的影響。

1 心臟運(yùn)動(dòng)模型

心臟運(yùn)動(dòng)為近似周期并且可以分解為心跳運(yùn)動(dòng)及呼吸運(yùn)動(dòng)2個(gè)主要頻率組成部分。基于心臟運(yùn)動(dòng)的以上特點(diǎn)，擬采用頻率域諧波合成方法進(jìn)行建模。諧波運(yùn)動(dòng)模型是以心臟數(shù)據(jù)頻域特征為基礎(chǔ)對(duì)心臟運(yùn)動(dòng)的一種近似。從傅里葉譜分析中可以觀測(cè)[2]，心臟運(yùn)動(dòng)功率譜由很多諧波構(gòu)成?？紤]算法的復(fù)雜度，這里將模型的階數(shù)選取為7，由此模型可以表示為:

式中:常數(shù)項(xiàng)c0通過心臟測(cè)量數(shù)據(jù)均值得到；am、bm為傅里葉級(jí)數(shù)系數(shù)。在任意時(shí)刻t，真實(shí)心臟運(yùn)動(dòng)與諧波模型運(yùn)動(dòng)變化方程則可以表示為:

式中:yt代表了3D心臟運(yùn)動(dòng)的真實(shí)值；ut表示3D諧波運(yùn)動(dòng)；△yt表示心臟運(yùn)動(dòng)值的增量；△ut表示諧波運(yùn)動(dòng)的增量。用ut近似表示yt，△ut來近似表示△yt，可以進(jìn)一步假設(shè):

在式（4）近似的基礎(chǔ)上，可以得到諧波運(yùn)動(dòng)的狀態(tài)方程:

式中:xt與xt-1分別為t與t-1時(shí)刻的心臟表面水晶振子的位置；εn為過程噪聲引起的不確定度，用均值為零的多維高斯分布表示。

2 心臟運(yùn)動(dòng)測(cè)量模型

圖1為超聲傳感器水晶振子分布圖，可以將得到的10個(gè)測(cè)量值描述為:

這些測(cè)量值將通過下式校準(zhǔn):

式中:xt表示參考坐標(biāo)系中心臟表面水晶振子的坐標(biāo)；qi（i=1，2，…，5)表示參考坐標(biāo)系中底座上水晶振子的坐標(biāo)。式中描述采用歐氏距離，展開如下:

圖1 超聲傳感器水晶振子幾何分布圖

式中:αt來自于超聲波傳感器的測(cè)量；是計(jì)算值。測(cè)量概率可以表示為。從偏移值歸一化直方圖可知，測(cè)量與計(jì)算之間的誤差基本符合正態(tài)分布。測(cè)量誤差的主要來源是晶體的幾何位置分布，因?yàn)橥ㄟ^超聲波測(cè)得的距離與2個(gè)水晶振子中心的真實(shí)距離存在一定偏差。首先，水晶振子是用醫(yī)用膠粘在塑料底盤上的，因此并不能保證每個(gè)水晶振子的角度一致，從而由于角度的些微偏差而造成的幾何分布誤差很難測(cè)量。其次，水晶振子的直徑會(huì)給超聲測(cè)距帶來偏置，其誤差范圍在-0.5 mm～2 mm，最大的誤差出現(xiàn)在0 mm稍偏右處。基于這個(gè)事實(shí)，可以通過高斯分布對(duì)噪聲建模，模型的均值μ和方差σ的選取由范圍內(nèi)誤差值頻率決定。其余部分為隨機(jī)噪聲，可以解釋為測(cè)量系統(tǒng)中超聲回波產(chǎn)生的噪聲，可以采用均勻分布對(duì)其進(jìn)行建模。在整個(gè)誤差范圍內(nèi)分布。將這2部分分布疊加，并對(duì)每部分配權(quán)重系數(shù)whit和wrand，whit+wrand=1。測(cè)量概率可以進(jìn)一步被表達(dá)成以下形式:

那么，最終的概率分布則是2個(gè)分布的線性組合，如圖2所示。

圖2 超聲測(cè)量系統(tǒng)測(cè)量概率分布圖

3 貝葉斯濾波算法的兩種實(shí)現(xiàn)

3.1 擴(kuò)展卡爾曼濾波器

擴(kuò)展卡爾曼濾波器可以處理非線性模型問題[5]，因?yàn)槠浔硎镜氖菍?duì)真實(shí)信度的一種近似（一般通過高斯過程），也就是通過其均值和方差來近似表示其信度。在擴(kuò)展卡爾曼濾波中，狀態(tài)轉(zhuǎn)移概率和測(cè)量概率由非線性方程g、h分別表示[6]。

式中:εt和δt分別表示模型噪聲和測(cè)量噪聲。非線性方程g、h的線性化是由泰勒級(jí)數(shù)展開后取其線性項(xiàng)。

3.2 粒子濾波器

粒子濾波器是實(shí)現(xiàn)貝葉斯濾波的另一種方法[,7]。在擴(kuò)展卡爾曼濾波算法中，其通過高斯分布來近似表示狀態(tài)的分布；而在粒子濾波算法中，它是通過有限數(shù)目的采樣值（也稱為粒子）來近似表示后驗(yàn)概率分布。采樣值代表了后驗(yàn)分布，其可以表示為:

其中，M為粒子在粒子集合Xt中的總個(gè)數(shù)。任一狀態(tài)，代表了單一假設(shè)檢驗(yàn)下可能的真實(shí)狀態(tài)。粒子濾波器的核心思想是通過粒子集中有限數(shù)量粒子的組合來表示信度分布。

4 心臟跳動(dòng)數(shù)據(jù)實(shí)時(shí)濾波及性能驗(yàn)證

研究以提高整個(gè)手術(shù)輔助機(jī)器人系統(tǒng)的跟蹤精度為目的，提出采用狀態(tài)估計(jì)方法對(duì)心臟運(yùn)動(dòng)測(cè)量數(shù)據(jù)進(jìn)行實(shí)時(shí)濾波處理。在性能驗(yàn)證中這里用離線濾波處理數(shù)據(jù)做為標(biāo)桿數(shù)據(jù)，用實(shí)時(shí)數(shù)據(jù)對(duì)比離線數(shù)據(jù)作出算法性能結(jié)果。其中，粒子濾波算法中使用500個(gè)粒子來近似概率分布。

實(shí)驗(yàn)中將擴(kuò)展卡爾曼濾波算法（EKF）和粒子濾波算法（PF）分別應(yīng)用到布朗運(yùn)動(dòng)模型和諧波運(yùn)動(dòng)模型，布朗運(yùn)動(dòng)是對(duì)隨機(jī)過程的理想化近似，也稱為維納隨機(jī)過程，這里用它來表示跳動(dòng)心臟的運(yùn)動(dòng)。布朗模型更詳細(xì)的定義參見文獻(xiàn)[8]。圖3展示了基于諧波模型的擴(kuò)展卡爾曼濾波算法與粒子濾波算法得出的在線濾波數(shù)據(jù)和原始數(shù)據(jù)。實(shí)時(shí)濾波算法實(shí)驗(yàn)結(jié)果如表1所示。

圖3 3D心臟測(cè)量數(shù)據(jù)EKF與PF濾波處理前后波形比較圖

表1 實(shí)時(shí)濾波算法實(shí)驗(yàn)結(jié)果

從表1可看出，擴(kuò)展卡爾曼濾波算法和粒子濾波算法均可以有效地對(duì)心臟跳動(dòng)測(cè)量數(shù)據(jù)進(jìn)行實(shí)時(shí)濾波。從模型角度而言，諧波模型優(yōu)于布朗運(yùn)動(dòng)模型。從算法計(jì)算復(fù)雜度而言，擴(kuò)展卡爾曼濾波算法計(jì)算60 s心臟數(shù)據(jù)的時(shí)間要遠(yuǎn)小于粒子濾波算法。綜合而言，研究認(rèn)為擴(kuò)展卡爾曼濾波是最適合應(yīng)用到基于概率統(tǒng)計(jì)方法的實(shí)時(shí)濾波算法，因?yàn)槠錅?zhǔn)確性以及算法復(fù)雜度均適用于實(shí)時(shí)的手術(shù)輔助機(jī)器人系統(tǒng)。

5 結(jié)束語

在CABG手術(shù)輔助機(jī)器人系統(tǒng)中，高信噪比的測(cè)量是精確跟蹤控制的前提與關(guān)鍵。本文提出基于概率統(tǒng)計(jì)的實(shí)時(shí)濾波方法，首先建立了心臟運(yùn)動(dòng)的模型與測(cè)量系統(tǒng)的模型，并且分別得出帶有統(tǒng)計(jì)信息的狀態(tài)轉(zhuǎn)移概率與測(cè)量概率，應(yīng)用基于遞推狀態(tài)估計(jì)貝葉斯算法的擴(kuò)展卡爾曼濾波算法和粒子濾波算法實(shí)現(xiàn)了實(shí)時(shí)濾波方法，實(shí)驗(yàn)結(jié)果表明，基于概率統(tǒng)計(jì)的在線濾波方法不但優(yōu)于傳統(tǒng)的濾波方法，而且可以達(dá)到了傳統(tǒng)頻域?yàn)V波無法達(dá)到的濾波效果，將測(cè)量誤差控制在1.54 mm內(nèi)，為實(shí)際的手術(shù)輔助機(jī)器人控制算法在手術(shù)現(xiàn)場(chǎng)試驗(yàn)提供了可能。

［1］CAVUSOGLU M C，ROTELLA J，NEWEMAW S N，et al. Control algorithms for active relative motion cancelling for robotic assisted off-pump coronary artery bypass graft surgery［C］//12th International Conference on ICAR’05.USA:IEEE，2005:431-436.

［2］梁帆，孟曉風(fēng)，董登峰.心臟動(dòng)脈旁路手術(shù)中手術(shù)輔助機(jī)器人的模型跟隨控制［J］.光學(xué)精密工程，2012，20（1）: 131-139.

［3］BEBEK O，CAVUSOGLU M C.Intelligent control algorithms forroboticassistedbeatingheartsurgery［J］.IEEETransRobot，2007，23（3）:468-480.

［4］THRUNS，BURGARDW，F(xiàn)OXD.ProbabilisticRobotics［M］. New York:The MIT Press，2006.

［5］GRIMMETT G R，STIRZATER D R.Probability and Random Processes［M］.Oxford:Oxford University Press，2001.

［6］潘紹松，左洪福.工業(yè)增強(qiáng)現(xiàn)實(shí)中的相機(jī)跟蹤［J］.光學(xué)精密工程，2011，19（6）:1353-1359.

［7］龔俊亮，何昕，魏仲慧，等.采用改進(jìn)輔助粒子濾波的紅外多目標(biāo)跟蹤［J］.光學(xué)精密工程，2012，20（2）:413-420.

［8］GUBNER John A.Probability and Random Processes for Electrical and Computer Engineers［M］.Oxford:Cambridge University Press，2006.

Real-time filtering of beating heart motion based on Bayesian method

LIANG Fan，YANG Li-li，CUI Shi-gang，WU Xing-li
（Tianjin Key Laboratory of Information Sensing and Intelligent Control，Tianjin University of Technology and Education，Tianjin 300222，China）

A real-time filtering algorithm of cardiac motion measurement signal based on Bayesian filtering recursive state estimation is proposed in this paper to improve the accuracy of measurement of cardiac data，in which the error model of sensor and the model cardiac motion are established，and the realization of the method for extended Kalman filter and two Bayesian filtering methods for particle filter are adopted.The experiment results indicate that the proposed method greatly helps to enhance the quality of the measurement of the beating heart motion and provides the possibility of lowering the tracking error in the whole system.

surgical assisted robotics；real-time filtering method；extended Kalman filter；particle filter

R541；TN911.7

A

2095－0926（2015）02－0007－04

2014-12-08

國(guó)家自然科學(xué)基金資助項(xiàng)目（61178048）；天津市自然科學(xué)基金資助項(xiàng)目（14JCQNJC04300）；天津職業(yè)技術(shù)師范大學(xué)科研發(fā)展基金資助項(xiàng)目（KJY11-10，KYQD13022）.

梁 帆（1981—），男，講師，博士，研究方向?yàn)獒t(yī)療信號(hào)處理與手術(shù)機(jī)器人跟蹤控制、現(xiàn)代農(nóng)業(yè)智能環(huán)境控制以及智能控制等.