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

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 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 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 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.

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].

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

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. q₁ is selected as the origin of the coordinate frame. q₃ and q₅ form the x axis and y axis with the origin. d_ij 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:

Figure 4.

Geometrical configuration of the measurement system [22]

\begin{array}{l} d_{k 1}^{2} - {\hat{x}}_{t}^{2} = d_{k 3}^{2} - {(d_{13} - {\hat{x}}_{t})}^{2} \\ d_{k 1}^{2} - {\hat{y}}_{t}^{2} = d_{k 5}^{2} - {(d_{15} - {\hat{y}}_{t})}^{2} \\ {\hat{z}}_{t} = + \sqrt{d_{k 1}^{2} - {\hat{x}}_{t}^{2} - {\hat{y}}_{t}^{2}} \end{array}

(1)

Sample data of time domain ECG and heart motion signals in three axes is shown in Figure 5.

Figure 5.

ECG and heart motion signals in time domain

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.

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 ε

f (t) = \sum_{j = 1}^{p} A_{11} (j) f (t - j) + ε (t)

(2)

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 ε₁:

\begin{array}{l} f (t) = \sum_{j = 1}^{p} A_{11} (j) f (t - j) + \\ \sum_{j = 1}^{p} A_{12} (j) g (t - j) + ε_{1} (t) \end{array}

(3)

If the variance ε₁ 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]:

G C I_{g \to f} = \log \frac{ε_{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.

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. d[n], d̂[n] and e[n] represent the desired response, the 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.

Figure 8.

MVAR model based prediction algorithm diagram

3.1 MVAR model

The MVAR model mathematically represents Granger Causality. x₁(k) and x₂(k) are two sets of time series. We suppose that the temporal dynamics of x₁(k) and x₂(k) are represented by the bi-variate autoregressive processes as follows.

\begin{array}{l} x_{1} (k) = \sum_{j = 1}^{p} a_{11} (j) x_{1} (k - j) + \sum_{j = 1}^{p} a_{12} (j) x_{2} (k - j) + v_{1} (k) \\ x_{2} (k) = \sum_{j = 1}^{p} a_{21} (j) x_{1} (k - j) + \sum_{j = 1}^{p} a_{22} (j) x_{2} (k - j) + v_{2} (k) \end{array}

(5)

where p denotes the model order, a₁₂ and a₂₁ represent the cross-correlation coefficients, a₁₁ and a₂₂ represent the coefficients of auto-regression, andv1 and v₂ are white noise. If we combine the two time series, we obtain the MVAR model:

\begin{array}{l} \vec{x} (k) = A_{1} \vec{x} (k - 1) + A_{2} \vec{x} (k - 2) + \dots \\ + A_{p} \vec{x} (k - p) + \vec{v} (k) \end{array}

(6)

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

\begin{array}{l} x [k + 1] = F x [k] + B v [k] \\ y [k] = C x [k] + n [k] \end{array}

(7)

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

F = [\begin{matrix} A_{1} & A_{2} & \dots & A_{p} \\ I_{2} & 0 & \dots & 0 \\ 0 & I_{2} & ​ & ​ \\ ​ & ​ & ⋱ & 0 \end{matrix}]

(8)

State vector is defined as:

\begin{array}{l} x [k - 1] = [x_{1} (k - 1), x_{2} (k - 1), x_{1} (k - 2), \\ \dots, x_{1} (k - p), x_{2} (k - p)]^{T} \end{array}

(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.

3.2 Parameterization of MVAR model

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.

Table 1.

Equation formulation in dual KF [28]

Initialize with:

\begin{array}{l} {\hat{w}}_{0} = E [w] P_{w_{0}} = E [(w - {\hat{w}}_{0}) {(w - {\hat{w}}_{0})}^{T}] \\ {\hat{x}}_{0} = E [x]^{'} P_{x_{0}} = E [(x - {\hat{x}}_{0}) {(x - {\hat{x}}_{0})}^{T}] \end{array}

the time-update equations for the weight filter are

\begin{array}{l} {\hat{w}}_{k}^{-} = {\hat{w}}_{k - 1} \\ P_{w_{k}}^{-} = P_{w_{k - 1}} + R_{k - 1}^{r} = λ^{- 1} P_{w_{k - 1}} \end{array}

and those for the state filter are

\begin{array}{l} {\hat{x}}_{k}^{-} = F {\hat{x}}_{k - 1} \\ P_{x_{k}}^{-} = A_{k - 1} P_{x_{k - 1}} A_{k - 1}^{T} + R^{v} \end{array}

The measurement-update equations for the state filter are

\begin{array}{l} K_{k}^{x} = P_{x_{k}}^{-} C^{T} {({CP}_{x_{k}}^{-} C^{T} + R^{n})}^{- 1} \\ {\hat{x}}_{k} = {\hat{x}}_{k - 1}^{-} + K_{k}^{x} (y_{k} - {C \hat{x}}_{k}^{-}) \\ P_{x_{k}} = (I - K_{k}^{x} C) P_{x_{k}}^{-} \end{array}

and those for the weight filter are

\begin{array}{l} K_{k}^{w} = P_{w_{k}}^{-} {(C_{k}^{w})}^{T} {[C_{k}^{w} P_{w_{k}}^{-} {(C_{k}^{w})}^{T} + R^{e}]}^{- 1} \\ {\hat{w}}_{k} = {\hat{w}}_{k - 1}^{-} + K_{k}^{w} (y_{k} - {C \hat{x}}_{k}^{-}) \end{array}

where

C_{k}^{w} = C \frac{\partial {\hat{x}}_{k}^{-}}{\partial w} |_{w = \hat{w}}

Figure 9.

The dual extended Kalman filter. The top part estimates the states and the bottom part generates weight estimation.

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

\hat{d} [k] = \hat{F} [k] \hat{x} [k]

(9)

4. Experimental results and discussion

4.1 Experimental results

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.

Table 2.

Parameters comparison

Prediction Algorithm	Parameters	1st dataset	2nd dataset
AR-F	Filter order Forgetting factor	12 0.995	18 0.9999
FS-EKF	Number of harmonies Number of states	6 14	10 22
MVAR-F	Filter order of autocorrelation Filter order of cross-correlation Forgetting factor	12 12 0.995	12 18 0.9999

Table 3.

Prediction Results

Prediction RMS errors [μm]	AR-F	FS-EKF	MVAR-F
1st dataset
X axis	335	99.4	69
Y axis	334	56.1	36
Z axis	465	64.4	54
2nd dataset
X axis	223	120	102
Y axis	257	89.8	48
Z axis	375	63.7	33

Figure 10(a).

Heart motion prediction results in X axis using first dataset

Figure 10(b).

Heart motion prediction results in Y axis using first dataset

Figure 10(c).

Heart motion prediction results in Z axis using first dataset

Figure 11(a).

Heart motion prediction results in X axis using second dataset

Figure 11(b).

Heart motion prediction results in Y axis using second dataset

Figure 11(c).

Heart motion prediction results in Z axis using second dataset

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 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

\hat{y} (n) = \sum_{i = 1}^{M} w_{i} y (n - i) + v

(A1)

where the ŷ(n) is the estimation at time n of the heart motion, {w_i}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]:

\begin{array}{l} \hat{y} (k) = ψ {(k)}^{T} \hat{ω} (k - 1) \\ K (k) = \frac{P (k - 1) ψ (k)}{λ + ψ {(k)}^{T} P (k - 1) ψ (k)} \\ \hat{ω} (k) = ω (k - 1) + K (k) (y (k) - \hat{y} (k)) \\ P (k) = \frac{1}{λ} (P (k - 1) - \frac{P (k - 1) ψ (k) ψ {(k)}^{T} P (k - 1)}{λ + ψ {(k)}^{T} P (k - 1) ψ (k)}) \end{array}

(A2)

7. Appendix B: Time-varying Fourier series with Extended Kalman Filter

The time-varying Fourier series with a constant offset are given by:

y (t) = c + \sum_{i = 1}^{m} r_{i} (t) \sin θ_{i} (t)

(B.1)

where $θ_{i} (t) = i \int_{0}^{t} ω (τ) d τ + φ_{i} (t)$ , ω is the heart rate and $r_{i}$ and $φ_{i}$ are the amplitude and phase of the ith harmony. The state space model for this system is

\begin{array}{l} x (t + Δ t) = A (Δ t) x (t) + μ (t) \\ y (t) = h (x (t)) + v (t) \end{array}

(B.2)

with state vector defined as x(t) = [c(t),r_i(t),w(t),θ_i(t)]^T. Here the transfer matrix A could be derived if the variable c(t),r_i(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.

\begin{array}{l} P (t + Δ t | t) = A P (t | t) F^{T} + Q \\ S = σ_{R}^{2} + H P (t + Δ t | t) H^{T} \\ K = P (t + Δ t) H^{T} S^{- 1} \\ \hat{x} (t + Δ t | t + Δ t) = A \hat{x} (t | t) + K (z (t + Δ t) - h (A \hat{x} (t | t))) \\ P (t + Δ t | t + Δ t) = (I - K H) P (t + Δ t | t) \end{array}

(B.3)

A (Δ t) = [\begin{matrix} I_{m + 1} & ​ & ​ & ​ & 0 \\ ​ & 1 & ​ & ​ & ​ \\ 0 & Δ t & ⋱ & ​ & ​ \\ ​ & ⋮ & ​ & 1 & ​ \\ ​ & m Δ t & ​ & ​ & 1 \end{matrix}]

(B.4)

\begin{array}{l} H^{T} = {(\frac{δ h}{δ x})}^{T} | ​_{\hat{x} (t + Δ t) = A \hat{x} (t)} \\ = [\begin{matrix} 1 \\ \sin {\hat{θ}}_{1} (t + Δ t) \\ ⋮ \\ \sin {\hat{θ}}_{1} (t + Δ t) \\ 0 \\ {\hat{r}}_{1} (t + Δ t | t) \cos {\hat{θ}}_{1} (t + Δ t) \\ ⋮ \\ {\hat{r}}_{m} (t + Δ t | t) \cos {\hat{θ}}_{m} (t + Δ t) \end{matrix}] . \end{array}

(B.5)

Footnotes

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).

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.

References

Newman

M. F.

Kirchner

J. L.

Phillips-Bute

Gaver

Grocott

Jones

R. H.

. (2001) Longitudinal assessment of neurocognitive function after coronary artery bypass surgery. New England Journal of Medicine, vol. 344, no. 6, pp. 395–402.

Puskas

John D

Wright

Carolyn E

Ronson

Russell S

Brown

Morris W.

III Gott

John Parker

Guyton

Robert A

(1998) Off-pump multivessel coronary bypass via sternotomy is safe and effective. Ann. Thorac. Surg., vol. 66, pp. 1068–1072.

Trejos

A. L.

Salcudean

S. E.

Sassani

Lichtenstein

(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.

Cavusoglu

M.C.

Rotella

Newman

W. S.

Choi

Ustin

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.

Rotella

(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

Franke

T. J.

Bebek

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.

Franke

T. J.

Bebek

Cavusoglu

M. C.

(2008) Prediction of Heartbeat Motion with a Generalized Adaptive Filter, 2008 IEEE International Conference on Robotics and Automation Pasadena, CA, USA.

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

Thakral

Wallace

Tolmin

Seth

Thakor

(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.

10.

Ginhoux

Gangloff

J. A.

DeMathelin

M. F.

Soler

Leroy

Sanchez

M. M. A.

Marescaux

(2005) Active filtering of physiological motion in robotized surgery using predictive control. IEEE Trans. Robotics, vol. 21, no. 1, pp. 67–79.

11.

Richa

Bó

A.P.L.

Poignet

, (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.

Yuen

Shelten G.

Kettler

D. T.

Novotny

P. M.

Plowes

R R. D.

Howe

R. D.

(2009) Robotic Motion Compensation for Beating Heart Intracardiac Surgery. International Journal of Robotics Research, vol. 28, p. 1355.

13.

Yuen

Shelten G.

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

Renaud

Cuvillon

Laroche

Forgione

Gangloff

(2009) Motion Prediction for Computer-Assisted Beating Heart Surgery. IEEE Transactions on Biomedical Engineering, vol. 56, no. 11,pp. 2551–2563

15.

Ortmaier

(2003) Motion compensation in minimally invasive robotic surgery. Doctoral Thesis, Technical University of Munich, Germany.

16.

Cuvillon

Gangloff

de Mathelin

Forgione

(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

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

Sastry

(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

(2008) Robotic-Assisted Beating Heart Surgery, Ph.D. Thesis, Case Western Reserve University [Online] Available: http://robotics.case.edu/publications.html#thesis

23.

Wiener

(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

Chen

Bressler

S. L.

(2006) Granger causality: Basic theory and application to neuroscience. In: Schelter

Winterhalder

Timmer

(Eds.), Handbook of Time Series Analysis Wienheim: Wiley, pp. 438–460.

26.

Geweke

(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

(2001) Kalman Filtering and Neural Networks, John Wiley and Sons, p. 304.

31.

Haykin

(2001) Adaptive Filter Theory, 4th ed. Upper Saddle River, New Jersey: Prentice Hall.