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) j1
p
f(t) A11 (j)f(t j)
j1
p
(3)
A12 (j)g(t j) 1(t) j1
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
j1
j1
p
p
j1
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
dˆ
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
ˆ ww
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 mt 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