Multivariate Autoregressive Model Based Heart Motion ... - InTechOpen

3 downloads 0 Views 838KB Size Report
Abstract A robotic tool can enable a surgeon to conduct off-pump coronary artery graft bypass surgery on a beating heart. The robotic tool actively alleviates the.
ARTICLE International Journal of Advanced Robotic Systems

Multivariate Autoregressive Model Based Heart Motion Prediction Approach for Beating Heart Surgery Regular Paper

Fan Liang1,*, Xiaofeng Meng1 and Yang Yu2 1 Science Technology on Inertia Laboratory, Beihang University, Beijing, China 2 Information Systems and Quantitative Sciences Rawls College of Business at Texas Tech University, Lubbock, TX, USA * Corresponding author E-mail: [email protected]

  Received 24 Jun 2012; Accepted 9 Jan 2013 DOI: 10.5772/55771 © 2013 Liang et al.; licensee InTech. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract  A  robotic  tool  can  enable  a  surgeon  to  conduct  off‐pump  coronary  artery  graft  bypass  surgery  on  a  beating  heart.  The  robotic  tool  actively  alleviates  the  relative motion between the point of interest (POI) on the  heart surface and the surgical tool and allows the surgeon  to  operate  as  if  the  heart  were  stationary.  Since  the  beating  heart‘s  motion  is  relatively  high‐band,  with  nonlinear  and  nonstationary  characteristics,  it  is  difficult  to follow. Thus, precise beating heart motion prediction is  necessary  for  the  tracking  control  procedure  during  the  surgery.  In  the  research  presented  here,  we  first  observe  that Electrocardiography (ECG) signal contains the causal  phase  information  on  heart  motion  and  non‐stationary  heart  rate  dynamic  variations.  Then,  we  investigate  the  relationship  between  ECG  signal  and  beating  heart  motion  using  Granger  Causality  Analysis,  which  describes  the  feasibility  of  the  improved  prediction  of  heart motion. Next, we propose a nonlinear time‐varying  multivariate  vector  autoregressive  (MVAR)  model  based  adaptive prediction method. In this model, the significant  correlation  between  ECG  and  heart  motion  enables  the  improvement of the prediction of sharp changes in heart  motion  and  the  approximation  of  the  motion  with  www.intechopen.com

sufficient  detail.  Dual  Kalman  Filters  (DKF)  estimate  the  states and parameters of the model, respectively. Last, we  evaluate  the  proposed  algorithm  through  comparative  experiments using the two sets of collected vivo data.    Keywords  Heart  Motion  Prediction,  Multivariate  Autoregressive  Model,  ECG  Signal,  Dual  Kalman  Filter  Beating Heart Surgery 

  1. Introduction  Coronary  Artery  Bypass  Graft  (CABG)  is  a  surgical  procedure performed to reduce the risk of coronary heart  disease.  Conventional  CABG  with  a  cardio‐pulmonary  bypass  machine  allows  surgeons  to  operate  on  a  motionless  heart;  however,  the  patient  suffers  cognitive  loss  for  longer  [1],  and  consequently  post‐surgery  complications are increased in terms of hospital time and  cost [2]. Compared to the conventional approach, the new  off‐pump  CABG  is  preferred  because  the  heart  keeps  beating  during  the  surgery,  which  would  decrease  the  resulting  post‐surgery  complications.  However,  manual  J AdvYu: Robotic Sy, 2013, Vol. 10, 129:2013 Fan Liang, Xiaofeng Meng andIntYang Multivariate Autoregressive Model Based Heart Motion Prediction Approach for Beating Heart Surgery

1

tracking of the quick and complex beating heart motion is  beyond  human  ability.  Robotic  technology  provides  an  effective  alternative  solution  by  actively  alleviating  the  relative  motion  between  the  beating  heart  and  robotic  tools  such  that  the  surgeon  can  operate  on  a  relatively  stationary  surgical  scenario.  Technically,  the  robot  tool  synchronizes  the  beating  heart  with  surgical  equipment  to  realize  the  motion  compensation.  In  the  tracking  process,  sufficient  accuracy  and  fast  tracking  requirements are the two challenges for this robotics tool.  On  the  one  hand,  in  off‐pump  CABG  surgery  surgeons  operate  on  blood  vessels  with  a  diameter  ranging  from  0.5‐2.0  mm.  Considering  the  safety  guarantee,  the  dynamic position error has to be limited to within 1% of  blood  vessel  diameter  or,  alternatively,  in  the  order  of  100‐250  m  in  terms  of  root  mean  square  (RMS)  [4].  On  the  other  hand,  the  main  characteristic  of  beating  heart  motion is the relatively high bandwidth. The frequency of  the motion can be up to 26 Hz, requiring the robot tool to  follow  and  predict  the  POI  in  a  very  short  period.  The  purpose  of  this  research  is  to  improve  the  tracking  performance  of  the  surgery  robotics  through  developing  the  heart  motion  prediction  method  for  the  receding  horizon  model  predictive  controller  [4].  The  detailed  control algorithm will not be discussed here due to space  restriction.  The  motion  prediction  problem  is  to  estimate  future heart motion based on past observations. In Figure  1  the  green  spots  representing  the  past  observations  produce  the  future  horizon  predictions  in  black  via  a  certain prediction algorithm. 

  Figure 1. Heart motion prediction schematic diagram 

In  order  to  emphasize  ECG  signal  in  the  heart  motion  prediction  method,  we  categorize  the  prediction  algorithms  into  two  families:  heart  motion  prediction  without  involving  ECG  signal,  and  heart  motion  prediction involving ECG signal.    The  heart  motion  prediction  without  ECG  signal  can  be  roughly  separated  into  time  domain  methods  [5‐8]  and  frequency domain methods [9‐14]. The main limitation of  these current methods is that they do not explicitly model  nonlinear  factors  such  as  variation  between  different  2

Int J Adv Robotic Sy, 2013, Vol. 10, 129:2013

heartbeat  periods.  This  deficiency  in  the  mathematical  modelling  of  the  beating  heart  motion  dynamics  would  cause  reduced  accuracy  both  in  prediction  and  tracking  control.  In  order  to  overcome  this  limitation,  some  research  introduces  ECG  signal  into  the  prediction  algorithm.    For  the  ECG‐involved  prediction  methods,  the  common  feature  is  the  utilization  of  the  correlation  between  ECG  and heart motion. Ortmaier [15] finds a strong correlation  between  ECG  signal  and  heart  motion  trajectory.  He  points out that this relationship implies these two signals  could be used as interchangeable inputs in the prediction  algorithm.  Cuvillon  et  al.  [16]  apply  the  correlation  between heart motion velocity and the ECG signal in the  linear parameter variant Finite Impulse Response model,  in which the impulse at each QRS occurrence of the ECG  is  the  model  input.  Bebeke  and  Cavusoglu  [17]  employ  ECG in the estimation of heart motion. The characteristic  R  point  in  the  ECG  wave  is  detected  by  wavelet  transform  to  provide  the  varying  heartbeat  cycle  such  that  the  heartbeat  period  can  be  adjusted  online.  This  algorithm could reduce the prediction errors compared to  the  constant  heartbeat  period.  Duindam  and  Sastry  [18]  extract  the  cardiac  phase  and  frequency  from  ECG  to  estimate  the  future  heart  surface  motion.  Although  ECG  signal  is  intensively  investigated  for  prediction,  the  studies  mentioned  above  focus  on  extracting  certain  features from ECG. There is still no general mathematical  model that has been proposed to process cross‐correlation  between  ECG  and  other  time  series.  We  expect  more  useful  dynamic  ECG  information  could  be  involved  in  heart  motion  prediction  under  the  certain  mathematical  model framework.     In  this  paper,  following  the  time  domain  prediction  stream and taking advantage  of ECG signal, we propose  an  adaptive  MVAR  model‐based  prediction  approach,  which  utilizes  and  models  the  correlation  between  ECG  signal  and  beating  heart  motion  to  improve  the  heart  motion  prediction.  We  investigate  the  property  of  ECG  signal  and  find  that  1)  ECG  causally  precedes  beating  heart motion, and 2) ECG shares the same characteristics  in motion. With these two results, this proposed method  could  better  predict  the  nonlinearity  of  heart  motion  dynamics.  From  the  perspective  of  regression,  the  proposed method extends the AutoRegression (AR) model‐ based methods [6,7] by adding the cross‐correlation term  to  acquire  information  from  other  sources.  In  terms  of  ECG  signal,  we  extend  the  method  in  [15‐18]  by  innovatively modelling the correlation information as the  interdependent  term  in  the  MVAR  model.  The  time‐ varying MVAR model, which explicitly embeds the ECG  information,  has  the  ability  to  adapt  to  changes  in  heart  rate  and  sharp  changes  in  the  morphology  of  fast  heart  motion.   www.intechopen.com

This  paper  is  set  out  as  follows:  Section  2  characterizes  the  property  of  ECG  and  the  relationship  between  ECG  and  heart  motion  by  Granger  Causality  Analysis,  respectively.  Next,  we  formulate  the  adaptive  multivariate autoregressive prediction algorithm. Finally,  the  proposed  prediction  algorithm  is  evaluated  through  comparative  experiments  using  collected  vivo  heart  motion data.  2. ECG signal property and Granger Causality   2.1 ECG signal property   ECG  is  a  transthoracic  (across  the  thorax  or  chest)  interpretation of the electrical activity of the heart over a  period of time [19]. ECG measures the rate and regularity  of  heartbeats  by  a  series  of  waves  and  complexes  that  have been labelled in alphabetical order: the P wave, the  QRS wave, the T wave, and U wave, as shown in Figure 2  [20]. 

2.2 ECG and heart motion signal demonstration  All the ECG and heart motion data are collected from an  animal model (adult calf. The ECG and heart motion data  are measured simultaneously, ECG signal being collected  from  the  analogue  data  acquisition  part  of  the  Sonomicrometry  system  [17].  The  sonomicrometer  measures  the  distances  within  the  soft  tissue  by  ultrasound  signals.  The  position  of  heart  motion  on  the  heart  surface  is  located  at  a  point  one  third  of  the  way  from  the  starting  point  of  the  Left  Anterior  Descending  (LAD) artery. 3D configuration of the heart motion can be  calculated  by  using  distances  among  the  crystals  in  the  measurement  system  [21].  Figure  4  shows  the  geometry  of the sensors in the ultrasound measurement system.  q1 is  selected  as  the  origin  of  the  coordinate  frame.  q 3 and  q 5 form the x axis and y axis with the origin.  dij denotes  the  distance  between  the  ith  crystal  and  jth  crystal,  and  d ki denotes the distance between the POI and ith crystal.  The position of the POI meets the following equations: 

d2k1  xˆ 2t  d 2k3  (d13  xˆ t )2  

d2k1  yˆ 2t  d 2k5  (d15  yˆ t )2                       (1)  2 zˆ t   d k1  xˆ t2  yˆ t2

  Figure 2. ECG waveform complex 

Bebeke[17] points out that by using ECG the performance  of  the  position  estimation  could  be  improved  because  wave  forms  in  ECG  are  the  results  of  physiological  processes  that  causally  precede  heart  motion.  The  time  relationship  between  action  potential  and  mechanical  force  developed  by  the  ventricular  muscle  is  shown  in  Figure 3. In physiology, force development in the muscle  follows  rapid  depolarization  of  a  cardiac  muscle  fibre  [19,20].  The  completion  of  repolarization  coincides  approximately  with  the  peak  force.  The  duration  of  contraction  parallels  the  duration  of  the  action  potential,  which is about 150 to 200 ms. This relationship is the root  of the correlation between ECG and beating heart motion.  We  conduct  further  exploration  of  the  correlation  using  Fourier Analysis and Granger Causality Analysis below. 

  Figure 3. Time relationship between action potential and  mechanical force  www.intechopen.com

  Figure 4. Geometrical configuration of the measurement system  [22] 

Sample  data  of  time  domain  ECG  and  heart  motion  signals in three axes is shown in Figure 5.     Since  ECG  signal  causally  precedes  heart  motion,  the  Fourier  Analysis  result  reveals  there  are  the  common  characteristics between the two kinds of signals. It is clear  from Figure 6 that the main combinatorial frequencies in  both  ECG  and  heart  motion  signals  are  mostly  the  same  below  5  Hz.  This  result  tells  us  that  all  the  modes  that  dominate  the  ECG  are  inherited  in  the  heart  motion,  meaning  that  the  dynamics  of  ECG  could  be  utilized  to  infer the dynamics of heart motion.  Fan Liang, Xiaofeng Meng and Yang Yu: Multivariate Autoregressive Model Based Heart Motion Prediction Approach for Beating Heart Surgery

3

ECG signal

0

X position[mm]

-5 20

Y position[mm]

-500 20

10 5 0 -5 20

Z position[mm]

ECG amplitude[mv]

time domain ECG and heart motion signal

500

10 5 0 -5 20

21

22

23

24

25

26

27

5

28

Heart Motion X

0 21

22

23

24

25

26

27

28

Heart Motion Y

21

22

23

24

25

26

27

28

Heart Motion Z

21

22

23

24 time[s]

25

26

27

28

 

Figure 5. ECG and heart motion signals in time domain  ECG signal amplitude spectrum

ECG

15

ECG FFT

10

5

0 0

1

2

1.5

Heart M otion

3

4

5

6

7

8

3D heart motion amplitude spectrum heart motion x FFT heart motion y FFT heart motion z FFT

1

0.5

0 0

1

2

3

4 time[s]

5

6

7

8

 

Figure 6. Amplitude spectrum of ECG and Heart motion signal  

2.3 Granger Causality Analysis 

The main objective of the Granger Causality analysis of 3D  heart  motion  data  and  ECG  signal  is  to  explore  the  correlation  between  the  multivariate  time  series  and  utilize  this  relationship  to  achieve  better  heart  motion  prediction.  Granger  Causality  is  the  concept  adopted  and  formulized  from Wiener [23] by Granger [24]. Granger defines causality  in terms of predictability of time series, that is, if series g(t)  contains  information  in  past  terms  that  helps  in  the  prediction of series f(t), then g(t) is said to cause f(t) [25]. The  standard  implementation  of  Granger  Causality  is  via  the  MVAR model, in which a set of time series are modelled as  weighted sums of past values. More specifically, if we try to  predict a value of f(t) using p previous values of the f series  only, we get a prediction error      

p

f(t)   A11 (j)f(t  j)   (t)                          (2)  j1

p

f(t)   A11 (j)f(t  j) 

 

j1

p

                             (3) 

 A12 (j)g(t  j)   1(t) j1

If the variance   1 is lower than the variance   , g causes f  in the sense of Granger Causality. The Granger Causality  measure  in  time  domain  is  called  the  Granger  Causality  Index (GCI) [26]:   

GCIg  f  log

1                                   (4)  

Here  we  calculate  GCI  by  using  the  Matlab  Causal  connectivity analysis toolbox [27]. The resulting direction  graph  in  Figure  7  shows  us  the  causality  relationship  among them. 

If we try to predict a value of f(t) using p previous values  of  the  series  f,  and  p  previous  values  of  g,  we  obtain  another prediction error  1 : 

4

Int J Adv Robotic Sy, 2013, Vol. 10, 129:2013

www.intechopen.com

are represented by the bi‐variate autoregressive processes  as follows.  p

p

j1

j1

p

p

j1

j 1

x1 (k)   a11 (j)x1 (k  j)   a12 (j)x 2 (k  j)  v1 (k)

 

   (5) 

x 2 (k)   a 21 (j)x1 (k  j)   a 22 (j)x 2 (k  j)  v 2 (k)

  Figure 7. ECG signal and heart motion data causal direction 

As indicated in Figure 7, the ECG causes heart motion in  all  three  axes,  which  matches  the  observation  made  in  Section 2 that ECG signal causally precedes heart motion.  The  current  information  in  the  ECG  signal  could  imply  what  will  happen  next  in  the  beating  heart  motion.  Therefore,  we  include  both  ECG  and  beating  heart  information in MVAR.   3. MVAR model based prediction   In  order  to  track  the  beating  heart  precisely,  the  robot  controller  needs  the  prediction  of  heart  motion.  Figure  8  shows  the  prediction  scheme  of  immediate  future  heart  motion by using an adaptive time‐varying MVAR model.  ˆ and e[n]   represent  the  desired  response,  the  d[n] , d[n] predicted  response  and  the  error  signal  at  the  current  time,  respectively.  The  prediction  problem  is  to  approximate  the  beating  heart  motion  by  the  weighted  past  observations  of  itself  and  the  weighted  past  observations  of  ECG  signal.  The  MVAR  model  includes  the coupling elements rather than multiple channel auto‐ regressions.  The  DEKF  parameterizes  the  MVAR  model  whenever  new  data  are  added.  We  separate  four  sets  of  time series into three pairs of bi‐variant time series. Each  pair of bi‐variants includes ECG and each individual axis  heart motion data. 

where p denotes the model order,  a12 and  a 21  represent  the cross‐correlation coefficients,  a11 and a 22 represent the  coefficients  of  auto‐regression,  and v1   and  v 2   are  white  noise.  If  we  combine  the  two  time  series,  we  obtain  the  MVAR model:   

a (j) a12 (j)   where x(k)  [x1(k),x 2 (k)]T , A j   11  , j  1, ,p , a 21(j) a 22 (j)   v(k)  [v1(k),v 2 (k)]T  

  The state space canonical form of the MVAR model could  be easily reformulated as:    

d [n]

dˆ [ n ]

e[n ]

  Figure 8. MVAR model based prediction algorithm diagram 

3.1 MVAR model 

The  MVAR  model  mathematically  represents  Granger  Causality.  x1 (k) and  x 2 (k) are two sets of time series. We  suppose  that  the  temporal  dynamics  of  x1 (k) and  x 2 (k)   www.intechopen.com

x[k  1]  Fx[k]  Bv[k]                         (7)  y[k]  Cx[k]  n[k]

F  is the canonical form of dimensions 2p by 2p given by:  

 

 A1 A 2  A p    I 0  0                        (8)  F   2  0 I2     0 

State vector is defined as:   

 x[ n  1]

   x(k)  A1x(k  1)  A 2 x(k  2)                  (6)    A px(k  p)  v(k)  

x[k  1]  [x1 (k  1),x 2 (k  1),x1(k  2),  ,x1(k  p),x 2 (k  p)]T

         (9) 

In  this  model  the  coefficient  matrices  weight  the  information of the past of the entire multivariate system.  We model the causal interactions between ECG and heart  motion  by  the  off‐diagonal  elements  of  A j (j  1, ,p)   matrices. We model the influence of the history of a heart  motion  process  on  the  predicted  value  by  the  diagonal  elements.    With the MVAR model in hand, we implement each axis  heart motion and ECG by this bivariate model. In order to  generate  the  predictions,  we  should  estimate  all  the  coefficient  matrices A j   in  F.  However,  what  we  have  in  equation  (7)  is  the  measurement,  and  the  conversional  least square based adaptive algorithm cannot estimate the  matrix F directly. Therefore, we utilize the DKF [28‐30] in  the parameterization process.  Fan Liang, Xiaofeng Meng and Yang Yu: Multivariate Autoregressive Model Based Heart Motion Prediction Approach for Beating Heart Surgery

5

3.2 Parameterization of MVAR model 

4. Experimental results and discussion 

As  equations  (7)  and  (8)  indicate,  the  clean  state  is  not  always  available  and  the  parameters  of  the  time‐varying  MVAR  model  need  to  be  determined  recursively  whenever a new data are added. The DKF estimates both  the  states  and  parameters.  In  this  task  two  KFs  run  concurrently.  At  every  time  step,  a  state  filter  estimates  the  state  using  the  current  model  weight  estimation,  while  the  weight  filter  estimates  the  weights  using  the  current  state  estimation.  Figure  9  shows  a  schematic  depiction  of  the  prediction  system.  Table  1  gives  the  complete equations of the DKF. 

4.1 Experimental results 

xˆk 1

xˆk

xˆk



yk

wˆ k

wˆ k

wˆ k 1

  Figure 9. The dual extended Kalman filter. The top part  estimates the states and the bottom part generates weight  estimation.  Initialize with: 

ˆ 0 )(w  w ˆ 0 )T ] ˆ 0  E[w] Pw0  E[(w  w w ,   xˆ 0  E[x] P  E[(x  xˆ )(x  xˆ )T ] x0

0

We  evaluate  the  proposed  algorithm  by  comparing  two  other  prediction  methods  using  two  sets  of  prerecorded  data.  The  prediction  algorithm  processes  three  pairs  of  bivariate data simultaneously.    We  compare  the  proposed  method  (MVAR‐F)  with  the  time  domain  method  adaptive  AutoRegression  Filter  (AR‐ F)  and  the  frequency  domain  method  time‐varying  Fourier  Series  filter  (FS‐EKF)  to  highlight  that  the  correlation among the different time series could improve  the  prediction,  which  is  also  the  idea  of  Granger  Causality.  In  AR‐F  each  set  of  time  series  is  modelled  using  the  autoregressive  information  without  any  cross‐ correlation  information.  The  Recursive  Least  Square  (RLS) method updates the parameters of the AR filter (see  Appendix A) In FS‐EKF, beating heart motion time series  is the summation of the truncated Fourier series without  any  cross‐correlation  information.  Extended  KF  updates  the parameters of the FS filter (see Appendix B).     Table  2  summarizes  the  parameters  of  three  algorithms.  Table  3  presents  the  one‐step  prediction  errors  in  the  sense of root mean square. Figures 10‐11 show the results  curve using two datasets.    Prediction  Algorithm  AR‐F 

0

the time‐update equations for the weight filter are 

ˆ k  w ˆ k 1 w Pw k

 Pw

k 1

 R rk 1

  Pw

FS‐EKF 

 

1

k 1

and those for the state filter are 

xˆ k  Fxˆ k 1 Px  A k 1Px ATk 1  R v k

MVAR‐F 

 

k 1

The measurement‐update equations for the state filter are 

K xk  Px CT (CPx CT  R n )1 k

 Cxˆ k )

 

Prediction RMS errors [  m ] X axis  Y axis  Z axis 

and those for the weight filter are 

K kw

Pw (C kw )T [C kw Pw (C kw )T k k



ˆk w

ˆ k 1 w

where 

Ckw  C

xˆ k w

 K kw (y k

ˆ ww

e 1

R ]

 Cxˆ k )

 

6

X axis  Y axis  Z axis 

 

Table 1. Equation formulation in dual KF [28] 

2nd  dataset  18  0.9999  10    22  12    18    0.9999 

     

ˆ ˆ ˆ  F[k]x[k]                                   (9)  d[k]

 

Int J Adv Robotic Sy, 2013, Vol. 10, 129:2013

AR‐F  FS‐EKF 

1st dataset  335  334  465  2nd dataset  223  257  375 

MVAR‐F 

99.4  56.1  64.4 

69  36  54 

120  89.8  63.7 

102  48  33 

Table 3. Prediction Results 

Once  we  obtain  the  estimation  of  the  parameters  in  (8),  one‐step prediction could be calculated via:   

1st  dataset  Filter order  12  Forgetting factor  0.995  Number of harmonies  6  Number of states     14  12  Filter  order  of    autocorrelation  12  Filter  order  of  cross‐   correlation  0.995  Forgetting factor 

Table 2. Parameters comparison 

k

xˆ k   K xk (y k Px  (I  K xk C)Px k k xˆ k 1

Parameters  

www.intechopen.com

4

x ref xAR F xFS EKF xMVAR F

3

0.5 prediction errors[mm]

2 heart motion prediction[mm]

xAR F err xFS EKF err xMVAR F err

1

1

0

-1

0

-0.5

-1 -2 -1.5

-3

-4 18.2

18.4

18.6 time[s]

18.8

-2 18.2

19

18.4

18.6 time[s]

18.8

19

 

Figure 10(a). Heart motion prediction results in X axis using first dataset  6

2.5

y ref yAR F yFS EKF yMVAR F

5

2

4

1.5

3 prediction errors[mm]

heart motion position[mm]

yAR F err yFS EKF err yMVAR F err

2 1 0 -1

1 0.5 0 -0.5

-2

-1

-3 -1.5 -4 18.2

18.4

18.6 time[s]

18.8

19

18.2

18.4

18.6 time[s]

18.8

19

 

Figure 10(b). Heart motion prediction results in Y axis using first dataset  2.5

1.5

z ref zAR F zFS EKF zMVAR F

2

1

1.5 1 prediction errors[mm]

heart motion prediction[mm]

zAR F err zFS EKF err zMVAR F err

0.5 0 -0.5 -1 -1.5

0.5

0

-0.5

-1

-2 -2.5 18.2

18.4

18.6 time[s]

18.8

-1.5 18.2

19

18.4

18.6 time[s]

18.8

19

 

Figure 10(c). Heart motion prediction results in Z axis using first dataset 

www.intechopen.com

Fan Liang, Xiaofeng Meng and Yang Yu: Multivariate Autoregressive Model Based Heart Motion Prediction Approach for Beating Heart Surgery

7

2

1

xAR F err xFS EKF err xMVAR F err

0.8 0.6 0.4

0

prediction errors[mm]

heart motion prediction[mm]

1

x ref xAR F xFS EKF xMVAR F

-1

-2

0.2 0 -0.2 -0.4 -0.6

-3

-0.8 -4 18.2

18.4

18.6 time[s]

18.8

-1 18.2

19

18.4

18.6 time[s]

18.8

19

 

Figure 11(a). Heart motion prediction results in X axis using second dataset  4

y ref yAR F yFS EKF yMVAR F

3

0.8

2

0.6 prediction errors[mm]

heart motion position[mm]

yAR F err yFS EKF err yMVAR F err

1

1

0

-1

0.4 0.2 0 -0.2 -0.4

-2

-0.6 -0.8

-3

-1 -4 18.2

18.4

18.6 time[s]

18.8

19

18.2

18.4

18.6 time[s]

18.8

19

 

Figure 11(b). Heart motion prediction results in Y axis using second dataset  4

1

z ref zAR F zFS EKF zMVAR F

3

zAR F err zFS EKF err zMVAR F err

0.8 0.6 0.4 prediction errors[mm]

heart motion prediction[mm]

2 1 0 -1

0.2 0 -0.2 -0.4

-2 -0.6 -3 -4 18.2

-0.8 18.4

18.6 time[s]

18.8

19

-1 18.2

18.4

18.6 time[s]

18.8

19

 

Figure 11(c). Heart motion prediction results in Z axis using second dataset 

8

Int J Adv Robotic Sy, 2013, Vol. 10, 129:2013

www.intechopen.com

4.2 Discussion 

We  find  that  the  MVAR‐F  algorithm  has  the  best  prediction  accuracy.  We  evaluate  the  performance  of  the  MVAR‐F  algorithm  through  three  different  measures:  correlation  representation,  heartbeat  frequency  variation  and prediction accuracy    a) Correlation representation    For  AR‐F  and  FS‐EKF,  coupling  with  ECG  signal  is  not  taken  into  account.  However,  the  MVAR‐F  extends  the  AR  filter  by  considering  the  influence  of  the  cross‐ correlation  coefficients  between  the  bivariate  series.  In  this way, the MVAR model is more than a combination of  multichannel  AR  models.  The  advantage  of  MVAR‐F  over  AR‐F  and  FS‐EKF  is  that  it  can  take  into  account  nonlinearity in the dynamics of the beating heart motion  by  correlation  with  ECG  signal.  MVAR‐F  is  nonlinear  in  some sense due to the integration of ECG.    b) Heartbeat frequency variation    All the three methods take into account the slight changes  in the heartbeat frequency. For AR‐F, the coefficients vary  to fit a new model whenever the heartbeat rate varies. FS‐ EKF  changes  its  parameters  to  adapt  the  heart  rate  variation. The MVAR‐F works in the same way as AR‐F;  the  cross‐correlation  part  introduces  the  heartbeat  frequency  information  from  ECG  signal  into  the  MVAR  model, which helps the MVAR‐F update its parameters to  fit  the  new  model  in  both  auto‐correlation  and  cross‐ correlation.  Therefore,  the  adaptive  MVAR‐F  is  more  capable  of  catching  the  variation  of  the  heartbeat  frequency than AR‐F and FS‐EKF.    c) Prediction accuracy    From  the  RMS  error  results  in  Table  3,  we  can  note  that  the  MVAR‐F  improves  the  prediction  results  in  comparison  with  AR‐F  and  FS‐EKF.  The  reason  for  this  improvement  lies  in  the  correlation  in  the  adaptive  MVAR  model,  which  accounts  for  the  correlation  between  the  heart  motion  and  ECG  signal.  In  addition,  the  cross‐correlation  describes  the  non‐stationary  and  nonlinear behaviour of the heart motion data through the  ECG signal, which contains the time lead dynamics of the  heart  motion.  MVAR‐F  utilizes  all  common  internal  modes  inferred  in  ECG  to  predict  the  motion  of  beating  heart.     The  proposed  algorithm  needs  to  be  validated  in  vivo  with  a  model  predictive  controller.  First,  pre‐process  before  MVAR‐F  should  be  conducted  to  filter  ECG  and  calculate  the  3D  heart  motion  data.  Second,  in  order  to  save  more  computational  time  for  pre‐process,  MVAR‐F  could  be  made  more  computationally  efficient  by  www.intechopen.com

replacing DKF with signal Kalman Filter via rearranging  the  states  and  parameters  into  one  state  vector  –  though  this could make predictions less precise. To conclude, we  expect  more  robust  and  informative  prediction  with  the  MVAR‐F.  5. Conclusion  In  this  paper,  we  use  Fourier  Analysis  and  the  Granger  Causality  method  to  analyse  the  causal  relationship  between  the  certain  axis  heart  motion  and  ECG  signal.  The directed resulting causality could help to improve the  prediction of heart motion for beating heart surgery. We  propose  a  novel  nonlinear  adaptive  MVAR  model  based  prediction  method.  The  parameters  of  the  MVAR  model  are  identified  recursively  by  DKF.  The  results  from  the  comparative  experiments  show  that  the  proposed  algorithm  gives  convincing  results  and  has  good  properties.  In  terms  of  amplitude  coupling  with  other  biological signal, MVAR‐F has the cross‐correlation in its  structure  to  model  the  coupling  effect  and  make  better  predictions through it. In relation to heart rate variation,  it will change the coefficients of both the auto‐regression  part  and  cross‐correlation  part  to  adaptively  follow  the  time‐varying statistics of the beating heart motion.    In  the  future  we  will  introduce  other  biological  time  series  signals,  such  as  diaphragm  signal  to  represent  the  motion  of  respiration,  into  the  MVAR  model.  We  will  implement  the  MVAR‐F  prediction  with  the  control  strategy  in  the  3D  test  bed.  More  precise  tracking  performance of assisted robotics in beating heart surgery  is expected.  6. Appendix A: AR filter  The  mathematical  representation  of  heart  motion  in  a  certain axis could be expressed as   M

yˆ (n)   wi y(n  i)  v                         (A1) 

 

i 1

ˆ   is  the  estimation  at  time  n  of  the  heart  where  the  y(n) motion,  {wi } is  the  coefficient  vector,  y(n  i) is  the  observation of the heart motion, and v is the white noise.  M  is  the  order  of  this  filter.  The  estimation  is  obtained  online  using  a  recursive  least  square  (RLS)  algorithm  formulated as [31]: 

yˆ (k)   (k)T ˆ (k  1) P(k  1) (k)    (k)T P(k  1) (k) ˆ ˆ (k)   (k  1)  K(k)(y(k)  y(k)) K(k) 

 

P(k) 

1



(P(k  1) 

P(k  1) (k) (k)T P(k  1)

   (k)T P(k  1) (k)

    (A2)  )

Fan Liang, Xiaofeng Meng and Yang Yu: Multivariate Autoregressive Model Based Heart Motion Prediction Approach for Beating Heart Surgery

9

7. Appendix B: Time‐varying Fourier series   with Extended Kalman Filter  The time‐varying Fourier series with a constant offset are  given by: 

9. References 

m

y(t)  c   ri (t)sin i (t)                      (B.1) 

 

i 1

t

where   i (t)  i  ( )d  i (t) ,  is  the  heart  rate  and ri 0 and  i are  the  amplitude  and  phase  of  the  ith  harmony.  The state space model for this system is  

 

x(t  t)  A( t)x(t)   (t) y(t)  h(x(t))  v(t)

                   (B.2) 

with  state  vector  defined  as  x(t)  [c(t),ri (t),w(t),i (t)]T .  Here the transfer matrix A could be derived if the variable  c(t),ri (t), w(t),i (t)  could evolve through a random walk.  Also  the  h(x(t)) could  be  formulated  as  y(t) in  the  measurement  equation  corrupted  by  measurement  uncertainty  variable  v(t) .  The  Extended  Kalman  Filter  (EKF) is used for state estimation. The detailed derivation  of the EKF for this model could be formulated as. 

P(t  t|t)  AP(t|t)FT  Q S   R2  HP(t  t|t)HT  

K  P(t  t)HT S 1 (B.3)  ˆx(t  t|t  t)  Ax(t|t) ˆ ˆ  K(z(t  t)  h(Ax(t|t))) P(t  t|t  t)  (I  KH)P(t  t|t)

 

I m  1 0   1                  (B.4)  A( t)   0 t     1    mt 1  

HT  (

 

Thanks  to  Dr  Cavosoglu  of  Case  Western  Reserve  University,  Cleveland,  US  for  help  and  advice  in  the  research relating to this paper during Fan Liang’s visiting  period in the US. 

h T ) |ˆ t)  Ax(t) ˆ  x x(t

  1   ˆ  sin (t t)     1      .                (B.5)    sin ˆ1 (t  t)   0    rˆ (t  t|t)cosˆ (t  t)  1 1        ˆ rˆm (t  t|t)cos m (t  t)   

8. Acknowledgments   This  study  was  partly  supported  by  the  China  National  Science Foundation for Group Creative Projects (number  61121003)  of  the  Science  Technology  on  Inertia  Laboratory (Beihang University, Beijing, China).  10 Int J Adv Robotic Sy, 2013, Vol. 10, 129:2013

[1]  Newman M. F., Kirchner J. L., Phillips‐Bute B., Gaver  V., Grocott H., Jones R. H., et al. (2001) Longitudinal  assessment of neurocognitive function after coronary  artery  bypass  surgery.  New  England  Journal  of  Medicine, vol. 344, no. 6, pp. 395–402.  [2]  John  D  Puskas,Carolyn  E  Wright.,Russell  S  Ronson.,W.Morris  Brown  III,  John  Parker  Gott,Robert  A  Guyton(1998)  Off‐pump  multivessel  coronary bypass via sternotomy is safe and effective.  Ann. Thorac. Surg., vol. 66, pp. 1068–1072.  [3]   Trejos  A.  L.,  Salcudean  S.  E.,  Sassani  F.,  Lichtenstein  S.  (1999)  On  the  feasibility  of  a  moving  support  for  surgery  on  the  beating  heart.  In:  Proc.  of  Medical  Image  Computing  and  Computer‐Assisted  Interventions  (MICCAI),  Cambridge,  UK,  pp.  1088– 1097.  [4]  Cavusoglu  M.C.,  Rotella  J.,  Newman  W.  S.,  Choi  S.,  Ustin  J.,  Sastry  S.  S.  (2005)  Control  algorithms  for  active  relative  motion  cancelling  for  robotic  assisted  offpump  coronary  artery  bypass  graft  surgery.In:  Proc.  of  the  12th  International  Conference  on  Advanced  Robotics  (ICAR),  Seattle,  WA,  USA,  July  2005, pp. 431–436.  [5]  Rotella  J.  (2004)  Predictive  tracking  of  quasi  periodic  signals  for  active  relative  motion  cancellation  in  robotic assisted coronary artery bypass graft surgery.  MS  Thesis,  Case  Western  Reserve  University,  Cleveland,  OH,  USA.  [Online].  Available:  http://robotics.case.edu/publications.html#thesis  [6]  Franke  T.  J.,  Bebek  O.,  Cavusoglu  M.  C.  (2007)  Improved  prediction  of  heart  motion  using  an  adaptive  filter  for  robot  assisted  beating  heart  surgery.  In:  Proc.  IEEE/RSJ  Int.  Conf.  Intell.  Robots  Syst., pp. 509–515.  [7]  Franke  T.  J.,  Bebek  O.,  Cavusoglu  M.  C.  (2008)   Prediction  of  Heartbeat  Motion  with  a  Generalized  Adaptive  Filter,  2008  IEEE  International  Conference  on Robotics and Automation Pasadena, CA, USA.  [8]  King  A.  P.,  Rhode  K.  S.,  Razavi  R.  S.,  Schaeffter  T.R.  (2009)  An  Adaptive  and  Predictive  Respiratory  Motion  Model  for  Image‐Guided  Interventions:  Theory  and  First  Clinical  Application.  IEEE  Transactions on Medical Imaging, vol. 28, no. 12, pp.  2020 ‐ 2032   [9]  Thakral A., Wallace J., Tolmin D., Seth N., Thakor N.  (2001)  Surgical  motion  adaptive  robotic  technology  (S.M.A.R.T):  Taking  the  motion  out  of  physiological  motion.  In:  Proc.  Int.  Conf.  Med.  Image  Comput.  Comput.‐Assisted  Intervention  (Lecturer  Notes  in  Computer Science, vol. 2208), pp. 317–325.  www.intechopen.com

[10] Ginhoux  R.,  Gangloff  J.  A.,  DeMathelin  M.  F.,  Soler  L.,  Leroy  J.,  Sanchez  M.  M.  A.,  Marescaux  J.  (2005)  Active  filtering  of  physiological  motion  in  robotized  surgery  using  predictive  control.    IEEE  Trans.  Robotics, vol. 21, no. 1, pp. 67–79.  [11] Richa,  R.,  Bó,  A.P.L.  ;  Poignet,  P.  ,  (2010)  Beating  Heart Motion Prediction for Robust Visual Tracking.  2010  IEEE  International  Conference  on  Robotics  and  Automation  Anchorage  Convention  District,  Anchorage, Alaska, USA.pp. 4579 ‐ 4584  [12] Shelten  G.  Yuen,  Kettler  D.  T.,  Novotny  P.  M.,  R  Plowes  R.  D.,  Howe  R.  D.  (2009)  Robotic  Motion  Compensation  for  Beating  Heart  Intracardiac  Surgery.  International  Journal  of  Robotics  Research,  vol. 28, p. 1355.  [13] Shelten  G.  Yuen,  Perrin  D.  P.,  Vasilyev  N.  V.  (2010)  Force  Tracking  with  Feed‐Forward  Motion  Estimation  for  Beating  Heart  Surgery.  IEEE  Transactions on Robotics, vol. 26, no. 5,pp. 888 ‐ 896  [14] Bachta  W.,  Renaud  P.,  Cuvillon  L.,  Laroche  E.,  Forgione A., Gangloff J. (2009) Motion Prediction for  Computer‐Assisted  Beating  Heart  Surgery.  IEEE  Transactions on Biomedical Engineering, vol. 56, no.  11,pp. 2551 ‐ 2563  [15] Ortmaier  T.  (2003)  Motion  compensation  in  minimally  invasive  robotic  surgery.  Doctoral  Thesis,  Technical University of Munich, Germany.  [16] Cuvillon L., Gangloff J., de Mathelin M., Forgione A.  (2005)  Toward  robotized  beating  heart  TECABG:  Assessment  of  the  heart  dynamics  using  high‐speed  vision.  In:  Proc.  Int.  Conf.  Med.  Image  Comput.  Comput.‐  Assisted  Intervention  (Lecturer  Notes  in  Computer Science, vol. 3750), pp. 551–558.  [17] Bebek O., Cavusoglu M. C. (2007) Intelligent control  algorithms for robotic‐assisted beating heart surgery.  IEEE Trans. Robotics, vol. 23, no. 3, pp. 468–480.  [18] Duindam  V.,  Sastry  S.  (2007)  Geometric  motion  estimation  and  control  for  robotic  assisted  beating‐ heart  surgery.  In:  Proc.  of  IEEE/RSJ  International  Conference  on  Intelligent  Robots  and  Systems  (IROS), San Diego, CA, USA, pp. 871–876.   [19] Berne  R.  M.,  Levy  M.  N.  (2001)  Cardiovascular  Physiology, 8th ed. St. Louis, MO: Mosby. 

[20] Berne  R.  M.,  Levy  M.  N.  (2000)  Principles  of  Physiology, 3rd ed. St. Louis, MO: Mosby.  [21] Ratcliffe M. B., Gupta K. B., Streicher J. T., Savage E.  B.,  Bogen  D.  K.,  Edmunds  J.  L.  H.  (1995)  Use  of  sonomicrometry  and  multidimensional  scaling  to  determine  the  three‐dimensional  coordinates  of  multiple  cardiac  locations:  Feasibility  and  initial  implementation.  IEEE  Trans.  Biomed.  Eng.,  vol.  42,  no. 6, pp. 587–598.  [22] Bebek  O.  (2008)  Robotic‐Assisted  Beating  Heart  Surgery,  Ph.D.  Thesis,  Case  Western  Reserve  University [Online] Available:     http://robotics.case.edu/publications.html#thesis  [23] Wiener  N.  (1956).  The  theory  of  prediction.  In:  Beckenbach  E.  F.  (Ed.)  Modern  Mathematics  for  Engineers, Chap 8. New York : McGraw‐Hill.  [24] Granger C. W. J. (1969). Investigating causal relations  by  econometric  models  and  cross‐spectral  methods.  Econometrica, vol. 37, pp. 424–438.  [25] Ding  M.,  Chen  Y.,  Bressler  S.  L.  (2006)  Granger  causality:  Basic  theory  and  application  to  neuroscience.  In:  Schelter  S.,  Winterhalder  M.  Timmer J. (Eds.), Handbook of Time Series Analysis  Wienheim: Wiley, pp. 438–460.  [26] Geweke J. (1982) Measurement of linear dependence  and  feedback  between  multiple  time  series.  J.  Am.  Stat. Assoc., vol. 77, no. 378, pp. 304–324.  [27] Seth  A.  K.  (2010).  A  MATLAB  toolbox  for  Granger  causal connectivity analysis. Journal of Neuroscience  Methods, vol. 186, pp. 262–273.  [28] Wan E. A., Nelson A. T. (1997) Neural dual extended  Kalman  filtering:  applications  in  speech  enhancement  and  monaural  blind  signal  separation.  In:  Neural  Networks  for  Signal  Processing  of  the  1997 IEEE Workshop, pp. 466–475.  [29] Wan  E.  A.,  Nelson  A.  T.  (2002)  Dual  Extended  Kalman Filter Methods, John Wiley & Sons, pp. 123‐ 173:  [30] Haykin  S.  (2001)  Kalman  Filtering  and  Neural  Networks, John Wiley and Sons, p. 304.  [31] Haykin  S.  (2001)  Adaptive  Filter  Theory,  4th  ed.  Upper Saddle River, New Jersey: Prentice Hall.   

 

www.intechopen.com

Fan Liang, Xiaofeng Meng and Yang Yu: Multivariate Autoregressive Model Based Heart Motion Prediction Approach for Beating Heart Surgery

11